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(i) PREFACE 


This volume is the third of a three volume set presenting 
the description and progi*am documentation of a mathematical 
model package for thermal pollution analyses and prediction. 

Two sets of programs, both in the free-surface formulation, 
are presented and clearly explained in this volume. These 
programs were developed by the Thermal Pollution Group at the 
University of Miami, and were funded by NASA, thus the program 
names NASUM II and NASUM III were given to reflect this joint 
effort. 

These models are three-dimensional and time dependent 
using the primitive equation approach. They have sufficient 
generalality in programing procedure to allow application at 
sites with diverse topographical features. Both programs pre- 
dict surface height variations i velocity field and temperature 
field for the "complete field". In the case of NASUM II a far- 
field formulation is used without including the plant thermal 
discharge; and in the case of NASUM III, a horizontal stretching 
is used to take account of the plant thermal discharge, and also 
to include far- field Influences such as varying tide and ambient 
currents at points sufficiently far from the point of discharge. 

These volumes are intended as user's manuals and, as such, 
present specific instructions regarding data preparation for 
program execution and specific simple problems. 
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(ii) LIST OF SYMBOLS 

The following list of symbols which are obtained from Volume 1 
are presented here for convenience. 

First term in n^(t), which is defined below 
A» Coefficient of second term in n (t) 

4 0 

Cq phase velocity of surface gravity waves, •'{gH* 

Cp specific heat at constant pressure 

f Coriolis parameter 

g acceleration due to gravity 

h depth relative to the mean water level 

H depth contour relative to free surface, h + n 

I grid index in x-direction or a direction 

J grid index in y-direction or 3 direction 

K grid index in z-direction or a direction 

k thermal conductivity 

K horizontal eddy viscosity 

H 

vertical eddy viscosity 
Kg surface heat transfer coefficient 

1 width of bay at ocean-bay interface 

L horizontal length scale 

P pressure 

Pg surface pressure 
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T Tempera cure 

T^ir Air temperature 

7 Water ambient temperature 

amb 

Equilibrium temperature 
t time 

t^ time lag in n^(t), which is defined below 

u velocity in x-direction (dimensional) 

V velocity in y-direction (dimensional) 

v^ Amplitude of inlet tidal velocity v^(t) 

w Velocity in z-direction (dimensional) 

X Horizontal coordinate 

y Horizontal coordinate 

z Vertical coordinate 

Greek Letters 

a horizontal coordinate in stretched system, 

3 horizontal coordinate in stretched system, 

a vertical coordinate in stretched system, ■ 

p density 

(j) phase angle in tidal current velocity 

n transformed vertical velocity 

u) angular frequency of tidal wave 

T. surface shear stress in x-direction 

KZ 

surface shear stress in y-direction 

X Z 


■ X 

- y 

z/H 
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n free surface elevation 

n^Ct) inlet tide level - Aj^ + A2 cos w(t+t^) 

Horizontal Stretching Parameters 

X horizontal stretching coordinate in x-direction 

Y horizontal stretching coordinate in /'direction 

X' dX 

3a 

X" A 

da2 

Y' « 

dS 

V" d^Y 

a the distance at which minimum step size is desired in 

x-direction (see transformation relation below) 

b the distance at which minimum step size is desired in 

/'direction (see transformation relation below) 

a,b,c^,C2>C2,c^, d and e are related and defined by the 
following relationships 

a ■ a + Cj^ Sinh' (c2(X-d)} 

8 - b + C3 Sinh' (c4(Y'e)} 
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1. INTRODUCTION 

This volume contains dascriptlons of the NaSuH II and 
NASUM III computer programs together with instruct' Ions on how tc 
operate these programs , As outlined in Volume I of this 
report, NASUM H is a phree-dimensional , time dependent, free* 
surface model intended for use in large domains whare rather 
coarse resolution is satisfactory. In NASUM. II the horizontal 
distance between nodes in the rectangular grid employed in the 
model is the same throughout the domain. NASUM III has the same 
basic three-dimensional, free-surface character as NASUM II but 
includes a form of "horizontal stretching" which provides fine 
resolution in some parts of th«r: domain and coarse resolution in 
other parts. The horizontal distance between nodes in the rec- 
tangular grid employed in the model is consequently a function 
of location in the domain. 

The program descriptions, associated algorithms, flow 
charts, program symbols, choice of input data, and sample 
problems for NASUM II and NASUMIII are contained herein for the 
ready access of the computer programs by the user. Note, that the 
governing equations, approximations, simplifying assximptions , 
and numerical methods of solution are presented in Voliome I. 

NASUMII (the far- field version of the free surface model) 
has been applied to South Biscayne Bay, Florida and is pre- 
sentee with computer results in Lee and Sengupta's (1977) report 
on Three-Dimensional Thermal Pollution Models and in a paper by 
Sengupta, Lee and Miller (1977). The South Biscayne Bay is a 
relatively shallow estuary with the principal driving mechanism 
being tidal flux at the ocean-bay interface, although wind 
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effaces are clearly evidenced in Che northern part of the bay. 

NASUM III (the horizontally stretched version of the free 
surface model) has been applied Co Hutchinson Island, St. Lucie, 
Florida, which Is a coastal site with a submerged d'^.scharge. 
Horlztonal stretching was used In order to obtain resolution 
In Che neighborhood of the discharge, while at the same time a 
large horizontal domain could be covered. If a constant grid 
size had been used, this would have required an excessively 
large number of grid points to cover the same extent of the 
boundaries of the domain. Therefore, in order, to circumvent 
this problem, a hyperbolic sine (SINK) stretching transformation 
was used in both the horizontally lateral and transverse directions , 
respectively, to obtain a small grid size in the neighborhood 
of the discharge and an Increasingly larger grid size as dis- 
tance from Che discharge point Increased. Waldrop and Farmer (1973) 
suggested a tangent (TAN) stretching transformation; however, 
the SINH stretching transformation was found to have advantages 
in this study. The details of Che comparison between the 
tangent and sinh stretching formulas are presented in Volume I. 

Lee and Sengupta (1977) and Tsai (1977) present the results of 
this Hutchinson Island Investigation, 

The effects of variable bottom topography, spatio-temporal 
free surface variations, surface heat transfer based on the 
equilibrium temperature concept introduced by Edinger and Geyer 
(1971), tide level variation, resultant ambient currents, and 
meteorological conditions have been factored into these models . 

In addition, turbulence has been modelled by using the eddy 
transport concept, and the effects of baroclinicity have been 
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included. Again, cha user should refer co Volume I for the 
complete mathematical formulations, approximations and assumptions, 
and the numerical methods of solution. 
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2. Program Descriptions 

This section presents the computer program algorithms and 
associated flow charts (in standard notation, c.f.. Murrill and 
Smith (1975) for the NASUM II and NASUM III. 

2 . 1 NASUM II (The Far Field M odel) 

2.1.1 Description of Program Algorithm 

The program algorithm for NASUM II is a follows: 

a) Integrate the surface height equation using forward- 
time, central-space differencing initially (FTCS) , and, there- 
after, central-time, central-space is used (CTCS) . 

b) Integrate the u-momentum equation using forward-time , 
central-space differencing initially (FTCS), and, thereafter, 
central-time , central-space is used (CTCS) with DuFort-Frankel 
differencing applied to the vertical momentum diffusion term as 
given, for example, in Roache (1972). 

c) Integrate the v-momentum equation using forward-time, 
central-space differencing initially (FTCS), and, thereafter, 
central- time , central-space is used (CTCS) with DuFort-Frankel 
differencing applied to the vertical momentum diffusion term. 

d) The equivalent vertical velocity, ri, is then computed 
by knowing H, u and v. The spatial integration is performed by 
applying Simpson's rule, (c . f . Crandell (1955)). 

e) The energy equation is then intregrated over time using 
forward- time , central-space (FTCS) throughout. 

f) The density is calculated from the equation of state. 

g) The pressure field is calculated from the hydrostatic 
equation using the trapezoidal rule for spatial integration. 

h) Then, the value of computed real time (or simulation 
time) is checked and the steps a) through g) repeated if so 
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desired. Reference to the flow chart presented in Fig. 1 will 
clarify this last step in the program algorithm. 

2.1.2 Flow Chart 

The flow chart for NASUM II is presented in Fig. 1. 

2 . 2 NASUM III (The horizontally stretched model ) 

2.2.1 Description of Program Algorithm 

The program algorithm for NASUM III is as follows : 

a) Integrate the surface height equation using forward- 
time, central-space differencing initially (FTCS) , and , 
thereafter, central-time, central-space is used (CTCS) . 

b) Integrate the u-momentum equation using forward-time, 
central-space differencing initially (FTCS), and, thereafter, 
central-time, central -space is used (CTCS) without DuFort - 
Frankel differencing applied to the vertical momentum diffusion 
'-hen, since the vertical diffusion term does not govern the 
time step value as it does for a shallow estuary. 

c) Integrate the v-momentum equation using forward- 
time, central-space differencing initially (FTCS), and, taere- 
afeer, central- time , central-space is used (CTCS). 

d) The equivalent vertical velocity, is then computed 
by knowing H, u and v. However, for a submerged discharge 

at the bottom of the basin is not zero (c.f. Volume I). The 
spatial integration is performed by applying the trapezoidal rule. 

e) The energy equation is then integrated over time using 
forward- time , central-space (FTCS) initially, and, thereafter, 
central- time , central-space (CTCS) is used. 

f) The density is calculated from the equation of state. 

g) The pressure field is calculated from the hydrostatic 
equation using the trapezoidal rule for spatial integration. 
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h) Thon, the value of oomputed i*eal time (or simulation 
time) is checked and the steps a) throup.h p,) repeated if so 
desired. Reference to the flow chart presented in I'ip. l: 
will clarify this last step in the propram alporithm. 

“ • I*’ low Chart 

The flow chart for NASUM III is presented in Fig. 2. 
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3. LIST OF PROGRAM SYMBOLS 

3.1 Symbols Common to Far-Field and Horizontally Stretched 
Model Programs 

This section presents in alphabetical order the program 
symbols, in FORTRAN language, and :ho:.r definition for those 
symbols which are common to the Far-Field and Horizontally 
Stretched Models. In many cases the definition is shortened 
by referring to algebraic symbols already defined in 
section (ii) of this volume. 

A 

A1 : first term in n^Ct) 

A2 : Coefficient of second term in n^Ct) 

B 

BH : By 
BV : B^ 

C 

Cl : Coefficient of inertia term in momentxam equations 

CC : Coefficient of Coriolis term in momentum equations 

CP : Coefficient of pressure term in momentum equations 

CH : Coefficient of horizontal diffusion term in momentum 

equations 

CV : Coefficient of vertical diffusion term in momentum 
equations 


D 
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D(I,J,K) 

: u(a,e,a) at t ■ t + 


DT 

: time step, At 


DX 

Grid size in a - direction, 

Aa 

DY 

: Grid size in @ * direction. 

A3 

DZ 

: Grid size in a - direction. 

Ao 

DUM 

: 0 ' 1 3 (Hu) + 3(Hv)> do 

ff ® 3^r^ 


DIHUTX 

3 (HuT) 
3 a 


DIHUUX 

: 3 (Huu) 

3a 


DIHUVX 

3 (Huv) 
3 a 


DIHUVY 

: 3 (Huv) 

33 


DIHUX 

: 3 (Hu) at k = k 

3 a 


DIHUXl 

: 3 (Hu) at k = k-1 

3 a 


DIHVTY 

; 3 (HvT) 

33 


DIHWY 

3 (Hw) 
33 


DIHVY 

: 3 (Hv) at k = k 

33 


DIHVYl 

; 3 (Hv) at k = k-1 

33 


DIPX 

: ^ 
3 a 


DIPY 

• 3P_ 

3 3 



DlTX 

D2TX 

DITY 

D2TY 

D2T2 

DlUX 

D2UX 

DlUY 

D2uy 

DlUWZ 

D2UZ 


3T 

5a 



3_T 

38 


3 T 



3 u 
3a 


3 u 
3a^ 

3u 


38 


2 

3 u 
2 

38 

3 (ufl) 
3 a 

2 

3 u 



3 V 
3 a 


DlVX 


2 

3 V 

3 01 ^ 
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D2VX 

DlVY 

D2VY 

DIVWZ 

D2VZ 

DIWTZ 


3 V 


33 

2 


3 V 

: 3 (v^) 


3 0 


2 



V 

T 


; 3 (J^T) 

3a 


Ed.J.F) : 

v(a, 3,cr) 

et t»t+At 

ETA(I,J) : 

n(a, 3) 


ETAX(I,J) : 

3_n 



3 a 


ETAYd, J) : 

3n 

33 


FF : f 



Gd.J.K) : 

u(a, 3, a) 

at t=*t 

GR : g 



H(I,J,K) : 
: 

u(a, 3) 
H(a,3) 

a t t=t 

a.t t=t in energy equation 

HDl)M(I,J) : 

- /^ { 3 (Hu)_ +.3 (Hv) } Ota 

0 3 a 3 8 


1 
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I 


i 


I 


HI(I.J): 
HN1(I,J) : 
HS: 
HT(I,J): 
HTD(I.J) : 
HTE(I,J) : 
HTMIN: 

HX(I,J) : 
HY(I,J): 

I: 

II: 

12 : 

13 : 

14 : 

15 : 

16 : 

IBAY: 


IHITE : 


IN: 

INLET: 


IRUN: 


h(o,S) 

H(a,0) at t»t + At In energy equation 

Ks 

H(a,0) at t"t - At in momentum equations 
H(a,d) at t*t in momentum equations 
H(a,0) at t*t + At in momentum equations 
Minimum value of H(a,0) in bay 


in tt - direction 
index for inlet along MAR-1* 
index for inlet along MAR-1 
index for outlet along MAR-2 
index for outlet along MAR-2 
for inlet along MAR-3 
for outlet along MAR-4 
Parameter which when specified either provides 
a constant time step for a shallow bay or for a 
deep bay. 

(-0 for shallow bay; -1 for deep bay) 

-0 for specifying initial surface 
-1 for not specifying initial surface 
Number of grid points in a - direction 
-1 for inlet along MAR-1 
-2 for inlet along MAR-3 
-o for first run 
-1 thereafter 


»h 
3 a 

3 h 

TT 


Index 

Lower 

Upper 

Lower 

Upper 

Index 

Index 


*NOTE: The MAR(I,J) matrix is explained in section 6.1.3. Fig. 3 

illustrates the location of II through 16 and JI through J6. 
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J J: Index in B - direction 

Jl: Index for inlet along MAR>I 

J2: Index for outlet along MAR>2 

J3: Lower index for inlet along MAR-3 

J4: Upper index for inlet along MAR-3 

J5 : Lower index for outlet along MAR-4 

J6: Upper index for outlet along MAR-4 

JN: Number of grid points in B * direction 

K K: Index in a- direction 

KH: Kjj 

KN: Number of grid points in a - direction 

KV: 

L L: Index of time cycle without energy equation 

LL: Index of time cycle with energy equation 

LN: Number of time cycles without energy equation 

LNI : Number of time cycles with energy equation 

M M: Parameter for either specifying (t) at the inlet or 

specifying Ct) at the inlet (-1 for for case) 

MAR(I,J): Numbering system for grid system - used to 

distinguish between different boundary finite difference 
schemes . 

P P(I,J,K): P(a.3,a) 

R RO(I,J,K): p(a,B,a) for variable density case 

RR: p(a,3,cj) for constant density case 

T T(I.J,K): T(a,3,a) at t“t 

TA: Te 
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U 

V 

w 


TAUXi 

TZX 

TAUY: 

Tzy 

TFLAT: 

Time interval from initial flat surface 
desired hour 

THT: 

Time interval from start*up to high tide 

TI: 

Initial temperature for isothermal bay 

TNI(I, J,K): 

T(a,0,a) at t«t + At 

TPH; 

Time lag for (t) 

TPHI: 

Time lag for (t) 

TT: 

TTOT + TTOTl 

TTOT: 

Total run time without energy equation 

TTOTl : 

Total run time with energy equation 

IKI.J.K): 

u(a,6,a) at t - t - At 

V(I,J,K): 

v(a,6,a) at t»t - At 

VO: 

amplitude of V^(t> at inlet 

W(I,J,K); 

0(a, a.a) 

WUD: 

1 ra(Hu) . 3(Hv)"\da 

V-t3 + 38 ■> 

WZ(I,J,K): 

o)(a, a, a) 


(n^-o) to Some 


3 . 2 Additional Symbols for Horizontally Stretched Model Program 
This section presents in alphabetical order the program symbols, 
in FORTRAN language, and their definition for those additional symbols 
for the Horizontally Stretched Model. Again, symbols not defined 
here have already been defined in section (ii) . 

A 



A: 

Value of 

(X 

- d)/c 

B 


B: 

Value of 

(Y 

- e)/C 

D 


DEEX 

: Value 

of 

<^1 


DEEY 

: Value 

of 

=3 
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U 


TT 


w 


X 


DELX; 

Grid size in 

0 - direction, Aa 

DE!V; 

Grid slz« in 

3 - direction, AS 

DHDX: 



DHDY: 

"3T ^ ^ ^ 


EEEX: 

Value of d 


EEEY: 

Value of e 



HK: K 


HTX: 

HTY: 


3H 

3H 

TF 


T(I,J,K): 

T(a, 3,o) at t*t - At 

TN(I,J,K): 

T(a,3,o) at t*t 

TF(I,J,K): 

T(a, 3,o) at t*t + At 

TAIR: 

Air temperature 

TAM: 

Ambient temperature of water body 


UM: u(ot,3,a) 

VM: v(a . 3 , <?) 

WH: w(a,3,a) 


XX- 


dX 
3H" 

XXX: A 


aot^ 
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4. MAIN PROGRAMS 

This ••ccion prettnte a d«c«il«d dcacrlpclon of cha main pro- 
grams for the NASUM II and NASUM III. Tha main programs thamsalvas 
appaar in Sactlon 7.1. 

4.1 NASUM II (Far-Fiald Main Program ) 

Tha following main program ouclina and associatad dascripcion 
is for tha far-fiald varsion of tha free surface modal, Tha main 
program nama is FMAIN, and appaars in Section 7.1. 

a) Specify nt'ibar of grid points, IN, JN and KN in PARAMETER 
statement (although tha geometry of the domain of solution under 
consideration will not cover all tha grid points; MAR(I,J)*0 covers 
range of grid points outside the domain of solution, where MAR(I,J) 
is constructed as shown in Fig. 5 for application to tha South 
Biscayna Bay) . 

b) Specify IRUN«0 or 1. The value 0 is used for the first run 
only, and 1 is used thereafter. 

For IRUN - 0, READ2 and INITIA arc used. 

For IRUN - 1, READl is used. 

c) Specify LN, LNl, M, INLET, IBAY, IHITE, II, 12, 13, 14, IS, 

16, Jl, J2, J3, J4, J5, J6, VO, TPH, TPHl, Al, A2, HTMIN, THT, 

TFLAT, Cl, CC, CP, CH. CV, GR, FF, RR, DX, DY, D2, KH, KV, BH. BV , 

TI. See section 3.1 for definition of these symbols, and refer to 
section 6.1.4, to follow, for a sample input of these parameters. 

d) Specify TAUX, TAUY, TA, and HS, as defined in section 3.1, 
each hour . 

e) Specify DT as defined in section 3.1. (in seconds) 

f) For L*l, TTOT-0.0: Energy equation is got coupled to the system 

of governing equations. The following subroutines are used: 

HEIGHT 

TIDE or VEL 

DATA 

UVVEL 

WVEL 

PRES 

ETT 

g) For L>1: Energy equation is not coupled to the system of gover- 

ning equations, but central-time is used now after the first time step 
has been executed (for L»l) . 
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HEILN 

TIDE or VEL 

DATA 

UWELN 

WVEL 

PRES 

ETT 

OLDHT 

OLDUV 

h) For LL>i, TT-TTOT + TTOTl > 0: Entrgy oquation is couplod to 

tht systtm of govtrnlng tquations. The following subroutlnee 

are i.sco: 

HEILN 

TIDE or VEL 

TIDAL 

DATA 

UVVELN 

WVEL 

PRES 

ETT 

OLDHT 

OLDUV 

TEMP 

OLDT 

1) After the final time cycle is computed, the following subroutines 

are used for printing and storing on magnetic tape: 

PKPARA 

PRETA 

PRUV 

WW 

PRW 

PRTEMP 

STORE 
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4.2 NASUM III (Horizontally Str<tched Main Program) 

The following main program outline and associated description 
is for the horizontally stretched version of the free surface model. 
The main program name is FMAIN, and appears in Section 7.1. 

a) Specify number of grid points, JN, JN and KN in PARAMETER 
statement for the *omain of interest. 

b) Read in all the data required and logic parameters IRUM, LN, 

Cl, CC, CP, CH, CV, GR, FF, RR, HK, DX, DY, DZ, KH, KV, BH, BV, TAUX, 
TAVY, TAIR, DT, DELX, DELY, DEEX, DEEY, EEEX, EEEY. See sections 3.1 
and 3.2 for definition of these symbols; and refer to Section 6.2.4, 
to follow, for a sample input of these parameters. 

c) Generate a two-dimensional matrix MAR(I,J) for locating the 
position of the points in the domain. 

d) Initialize all the necessary quantities, as defined in Section 
3.1; specify the discharge conditions and bottom topography, 

e) Convert the real vertical velocities W into the transformed 
sigma coordinate vertical velocities, fi. 

f) Calculate the horizontal stretching parameters, x| X",y',Y 
for the set of governing equations. 

g) Calculate the new predicted dependent variables . 

h) Store the data and new predicted dependent variables on magnetic 
tape and print out these values at the desired time step. 

i) For L»l, TTOT-0.0: Forward- time differencing is used. The 

following subroutines are used: 

HEIGHT 

UWEL 

TEMP 

PRES 

ETT 
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j) For L>1 or TTOT>DT: Central- time differencing is used. 

The following subroutines are used; 


HEILN 

UVVELN 

WVEL 

TEMPN 

PRES 

ETT 

OLDUVT 

OLDUV 


k) After the final time cycle is computed, the following subroutines, 
are used for printing and storing on magnetic tape: 

STORE 

PRPARA 

PRETA 

PRUV 

PRW 

PRTEM 


5 . INPUT DATA 


The data that is required for the execution of the main program in 
either NASUM II or NASUM III is called Input Data. The data required 
is listed in the order it appears in the respective programs, and the 
corresponding FORMAT (in FORTRAN language) is given corresponding to 
each data symbol. Section 5.1.1 lists the data input required for 
running NASUM II, and Section 5.2.1 lists the data input required for 
running NASUM III. The actual calculation required for several of 
the input data is given in Section 6.1.3 for NASUM II, and in 6.2.3 
for NASUM III. Note, the data input symbols have already been 
defined in Section 3 of this volume. 


5.1 NASUM iKFar-Field Model) 

The following number of computer data cards, with proper FORMAT, 
is now given ^ order as they appear in the main program in order to 
execute NASUM II (refer to Section 3.1 for definition of these FORTRAN 
symbols) . The data that must be calculated beforehand is given in 
Section 6.1.3. 
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5.1.1 DATA REQUIRED AND FORMAT 


CARD NO. 

DATA 

FORMAT 

1 

IRUN 

15 

2 

LN 

15 

3 

LNl 

15 

4 

M 

15 

5 

INLET 

15 

6 

IBAY 

15 

7 

IHITE 

15 

8 

11,12,13,14,15,16 

615 

9 

J1,J2,J3,J4,J5,J6 

615 

10 

V0,TPH,TPH1,A1,A2 

Free 

11 

HTMIN,TliT 

Free 

12 

CI,CC,CP,CH,CV 

Free 

13 

GR.FF.RR 

Free 

14 

DX,DY,DZ 

Free 

15 

KH,KV 

Free 

16 

BH,BV 

Free 

17 

TI 

Free 

18 

DTAUX(13)** 

Free 

19 

DTAUY(13) 

Free 

20 

DTA(13) 

Free 

21 

DHS(13) 

Free 

22 

DT 

Free 


**NOTE: 13 values of TAUX and TAUY and 13 values of TA and HS 

are read in for variation each hour. Note, the letter 
”D" proceeds previously defined symbols (Section 3.1), 
since this was necessary for computer convenience. 
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•5.2 NASUM III (Horizontally Stretched Model ) 

The following niomber of computer data cards, with proper 
FORMAT, is now given ^ order as they appear in the main program 
in order to execute NASUM III (refer to Section 3.1 and 3.2) for 
definition of these FORTRAN symbols) . The data that must be 
calculated beforehand is given in Section 6.2.3. 

5.2.1 DATA REQUIRED AND FORMAT 


CARD NO. 

DATA 

FORMAT 

1 

IRUN 

15 

2 

LN 

15 

3 

CI,CC,CP,CH,CV 

Free 

4 

GR,FF,RR,HK 

Free 

5 

DX,DY,DZ 

Free 

6 

KH,y.V,BH,BV 

Free 

7 

TAUX,TAUY 

Free 

8 

TAIR 

Free 

9 

DT 

Free 

10 

DELX 

Free 

11 

DELY 

Free 

12 

DEEX 

Free 

13 

DEEY 

Free 

14 

EEEX 

Free 

15 

EEEY 

Free 


26 


6. SAMPLE CASES 

The following sample cases will Illustrate and clarify to 
the user the proper choice of programs, subprograms (or subroutines), 
calculation of input parameters, sample input, and sample output for 
NASUM II and NASUM III. 

6.1 NASUM II (Far Field Model) 

6.1.1 Problem Statement - Application to Biscayne Bay 

Given Biscayne Bay, in Dade County, Florida, as an example 

application site, compute the surface heights, n, velocity field 
u, V, w, and the temperature distribution T at 2:00 P.M. knowing 
the meteorological data, the IR data base, and the tide data base 
for April 15, 1975. The IR data base is assumed synoptic at 2:00 P.M., 
the wind velocity and ambient temperature is known every hour, and 
the tide height, with respect to the mean water level, is known as 
a function of time at the ocean bay interface. 

Use the NASUM II far-field, free surface model program to 
obtain the desired results 

6.1.2 Choice of Subroutine Programs 

Section 4.1 is followed in order to choose the proper subroutine 
programs for this sample case. Sections 6.1,3 and 6.1.4 to follow next, 
will clearly illustrate what steps the user must follow in order to 
obtain the desired results for this sample case. 

6.1.3 Calculation of Input Parameters 

a) Construct a three-dimensional grid system for the Biscayne Bay. 

Fig. 4 illustrates the horizontal grid for the bay superimposed on 
the actual geometry of the domain of interest. The governing equations 
have been transformed into the a, 3 , a coordinate system, which maps 
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the variable depth basin into a constant depth basin. Then , 
depending on the desired resolution of vertical structure, the number 
of vertical grid points is selected, that is KN. The values of 
IN and JN are then selected with consideration of desired horizontal 
resolution versus computer storage and overall computation time. 

Then, the next step is to specify IN, JN, KN in the main program, 
as mentioned in section 4.1. For this application IN-34, JN-11, 

KN-5. 

b) Next, LN and LNl as defined in section 3.1 are specified. 

First, however, the time step DT is computed based on the criteria 
given in Volume I. For the Biscayne Bay, DT is computed by the 
vertical momentum diffusion criterion: 

DT - At < /K^ . 'm ^ C 60 96cm) ^ /Scm^/sec, 

Thus, DT - 10 sec. is chosen for this sample case, in order to 
ensure numerical stability, 

For starting the program at (t®o) = 0 for April 15, 1975 
LN- 1342, LNl-1380, since the program is started at &:26-a.m, and 
run without the energy equation to 10:10 am, at which time the 
IR data base is read in, as an initial condition, from sub- 
routine TIDAL, Then, from 10:10 am to 2:00 P.M. the energy equation 
is included. This procedure ignores the effect of density currents 
on the momentum and surface height equations from 6:26 am to 10; lOara. 
However, these density currents are quite small for the Biscayne 
Bay which is dominated by the wind and the tidal flux at the ocean- 
bay interface. 

c) The eddy viscosity coefficients have been estimated by applying 
the ”4/3 scaling law” to , previously known water basin values used 
by other researchers. 
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d) Next, Che matrix MAR(I,J) is constructed based on this 
particular grid system as shown in Fig . 5 . 

e) The depth matrix HI (I,J) is then constructed by specifying 

the depth below the mean water level, at each hcrizontaH 

grid point. 

f) The initial temperature matrix T (I,J,K) is constructed for 
the bay by first plotting the IR data base surface isotherms on 
the horizontal grid, and tnen interpolating to specify the 
temperature at each grid point (I,J), (See Fig. 6 ) The bay is 
shallow and well mixed vertically, hence, the vertical LemperaLure 
variation is initially set equal to zero. Subroutine TIDAL reads 
in from data cards T(I,J,K). 

g) The tidal current velocity ampli ude Vo is computed from the followin 
formula, asstuning a 90° phase shift between the tide height and 
tidal current velocity at the ocean-bay interface, (Ippen 1966): 


V » 1 - 
o V 

where 


2aCo \ i 2tt1 

-H — ) v-r- 

a - |n^l 

1 * 10A6 

A = C^T 
o 


» 37 cm for April 15, 1975 
= 4.3 X 10^ cm/ sec 
= 16 X lO^cm 


12C - 1.86 X 10' cm>> 1 

o 


<h> 190 cm, average at ocean-bay interface. 

Thus, Vo= 90 c"/sec 

h) TPH =3138.7 sec. is computed by letting Vo ( t) »0 at 10:10 A. M. (high 

tide where Vo(t) = Vo cos [(t + TPH)-^] and t=o at 8:00 A.M. on 4/15/75. 

TPHl = 7800 + FLAT = 13440 sec. (where 7800 sec = 8:00 A.M. to 10:00 A.M.) 

where TFLAT = 5640 sec., the time from Ho(t) = 0 at 6:26 A.M. for 
April 15, 1975 to the beginning of the Va(t) run at 8:00 AI'I on April 15, 
1975. 
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I) The equilibrium temperature, TA, and the equilibrium coefficient 
of surface heat transfer HS, are computed following Harleman and 
Stolzenbach (1973). Appendix C presents these formula. 

J) The wind stresses TAUX, TAUY are dc^.termined as shown in 
Appendix A. 


6.1.4 Sample Input 

The far- field solution is obtained by either specifying Vo(t) 
or nQ<t) at the ocean-bay interface. These two sample cases will 
now be given; 

6. 1.4.1 Velocity Case 


IRUN - 0 

LN ■ 360 (begin at 0800 EST, 4/15/75, and compute until 0900 EST) 
LNl-1 


M«1 (Vo(t) specified at inlet) 

INLET " 1 (ocean-bay interface along MAR*1) 

IBAY • 0 (shallow bay), i.e., constant time step of lOsec is used even 
in shallow regions - this was done to avoid inordinately small time 


steps during low water. 


IHITE <■ 1 (regression surface not read in at time of high tide) 


II, 12, 13, 14, 15, 16- 7, 16, 31, 33, 35, 35(refer to Fig. 3) 
Jl, J2, J3, J4, J5, J6- 11, 1, 12, 12, 12, 12(refer to Fig. 3) 


Vo, TPH, TPHl, Al, A2 --90, 3138.7, -13420, 11.8872, 37.1856 
HTMIN, THT - 60.96, 7800.0 
TFLAT - 0.0 


Cl, CC, CP, CH, CV-1., 1., 1., l.,l. 
GR, FF, RR - 980. , .00006, 1. 

DX, DY, DZ - 160,000., 160,000., 

KH, KV - 10,000. , 5. 

BH, BV - 10,000. , 5. 




.25 
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TI - 24.5 
DTAUX (1) - -.37 
DTAUY (1) - .15 
DTA (1) - 31.7 
DHS (1) - .00129 


DT - 10 

Next, IRUN - 1, LN - 420, LNl - 1 (lOilOaw) 

Next, IRUN - 1, LN - 1, LNl - 1380 (2:00pm) 

6. 1.4.2 Tide Height Case 
IRUN - 0 

LN » 264 (begin at 0626 est, 4/15/75 and compute unitl 0710 est) 
LNl- 1 

M » 2 (ngCt) specified at inlet) 

INLET * 1 (ocean-bay interface closing MAR-1) 

IBAY = 0 (shallow bay) 

IHITE » 1 (regression surface not read in at time of high tide) 
II, 12, 13, 14, 15, 16 (refer to Fig. 3) 

Jl, J2, J3, J4, J5, J6 (refer to Fig. 3) 

Vo, TPH, TPHl, Al, A2 - -90, 3138.7, -13420., 11.8872, 37.1856. 
HIMIN, THT - 60.96, 13420. 

TFLAT = 5640 


Cl, 

CC, 

CP, CH, CV 

= 1. . 

GR, 

FF. 

RR = 980. , 

.00006 
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DX, DY, DZ - 160,000.. 160,000., .25. 
KH, KV - 10, 000. , 5. 

BH, BV - 10,000, 5. 

TI - 24.5 
DTAUX(l) - -.37 
DTAUY(l) - .15 
DTA(l) - 31.7 
DHS (1) - .00129 


DT - 

10 


Next, 

IRUN - 1, 

LN 

Next, 

IRUN - 1. 

LN 

Next, 

IRUN - 1, 

LN 

Next, 

IRUN - 1, 

LN 


300, LNl - 1 C8 :00am; 

300, LNl - 1 (9; 00am) 

420, LNl - 1 (10:10am) 

1, LNl - 1380 (2:00pm) 
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6.1.5 Program Execution Procedure (NASUM II) 

This section describes Che procedure by which Che user 
executes the NASUM II program. 

a) Input Parameters ; The user must first follow the steps 
outlined in section 4.1, and become quite familiar with all the 
input parameters listed in section 5.1 and section 6.1.3. 

b) First Run ; In order to obtain surface heights and three- 
dimensional velocity and three-dimensional temperature, Che main 
program FMAIN is executed. In FMAIN there are two tape units. One 
is a READ unit designated as Unit 7. The other is a STORE unit des- 
ignated as unit 8. During the first run, ther is no need for unit 7, 
and unit 8 has to be provided to store results on a magnetic tape. 

c) Run Continuation ; For extending the results, the run has 

to be continued. The magnetic tape which was "unit 8" in the first 
run will now be read" Unit 7", for reading the previously stored 
results. Another magnetic tape is now to be provided as "unit 8" for 
storing the extended run results. The above procedure can be repeated 
until the results are obtained for the desired time. It is to be 
noted that for the first run IRUN*0, and for extended runs IRUN-1. 
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6.1.6 Sample Output 


Tha output from tha modal sample run is listed as follows t 

a) parameters 
b4 Surface heights 

c) Horizontal components of velocity 

d) Vertical velocity component 
•e) Temperatures 
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6.2. NASUM III (Horizontal Stretched Model) ; 

6.2.1. Problem Statement ; 

Given the Hutchinson Island, St Lucie, Florida as an 
example application site, compute the three dimensional 
velocity and temperature distribution for the following 
discharge and mateorologlcal conditions. 

Discharge volume from the 1 

Condensers of Power PlantJ* 363,000 6 p.m. 

Discharge Temperature « 35.0^C 

Ambient Temperature <■ 25.0^0 

Air Temperature - 30.0°C 

Current » 2 cm/ sec South 

6.3. Calculation of Input Parameters ; 

In this section, the specification of grid system, 
reference quantities and calculation of discharge velocities 
chosen will be presented first followed by the actual cal* 
culatlon of input data as they appear in the main program. 

6.3.1. Grid System 

The remote sensing data and ground truth data was 
available for the Hutchinson Island site and it is used 
to determine the size of the domain. The domain selected 
was 2380 m x 2000 m. The domain has a variable depth and 
so a variable bottom topography is used. A horizontally 
stretched grid system as shown in Fig.( ) is used. This 
would give more resolution of the plume in the near field 
where the effects of the thermal discharge are predominant. 




1 
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The grid system is also stretched in the vertical direction 
for ease of programming . In this way, same nxuiber of grid 
points can be used in shallow and deep regions of the basin. 
6.3.2. Calculation of Discharge Velocity ; 

In the numerical model a 9 point discharge is chosen. 
The discharge velocity is calculated by balancing the 

mass as shown below. 

For the numerical grid system, the mass into the domain * 
(Discharge Area x Velocity) ■ A x V 

The grid size chosen at the discharge is the minimum grid 
and equal to 50 m. So the discharge area is equal to 
(150 X 150) m^ as shown below 


1 : 

n 


1 

h— - - ( 

1 

‘ iO 

A 

^150m ^ 


A = (150 X 150) X 10“ cm^ 

V X A = Discharged volume (given) 

= 363,000 G.P.M 

V X (150x150x10“) = 363,000 x 0.0038 x 10 ^ cm^ 

515 sec 

V » 0 . 102 cm/ sec 

.‘.Velocity at discharge or inlet velocity for the model is 
0.102 cm/ sec. 

6.3.3. Reference Quantities : 

The reference eddy viscosity and diffusivity are de- 
termined using the 4/3 power law as below 

*re£ ■ 



•u 
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Where is the reference eddy vlecoclty and L it the 

maximum length of the domain in centimeters. 

Aref - 0.0025 (2380 x 
le 40.000 cm^/sec 

^ref 

For turbulent prandd number of 1 (Pr.. ■- « ■— -) 

^ ®ref 

®ref diffusivity) ■ 40,000 cm^ 

sec 

The vertical eddy viscosity and diffusivity chosen was 

10 cmf. 
sec 

6.3.4. Calculation of Incut Data as It Appears In 
the Main Program 

Card No Fortran Quantity 

— I — fmn — 

IRUN 0 for the first run and equal to 1 for later runs. 

Card No Fortran Quantity 

2 

LN is the number of cycles required. It is always advised 

to run the program for 10 or 13 cycles and check how the 

model is running. 

Card No Fortran Quantity 

3 kSTORfi 

If KSTORE Is equal to zero the model sill store the results 
on the tape to be provided and if it is eqtial to 1 the model 
will not store results on the tape. For KSTORE equal to 
zero or 1 the model will print results at the end of the 
rxin, 

Card No Fortran Quantity 

5 CI.CC.CH.CV.'CP 
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The values of these quantities are equal to 1,0 

Card No Fortran Quantity 

5 -TOTTROK 

GR Is gravity ■ 980,0 cm/ sec 

FF Is cor lolls term fi 0.0006 

For small domain this terms Is negligible and may be kept 
equal to 0.0 for all practical purposes 
RR Is density ■ 1.0 
HK - 

h - 1200 BTU/day °F - ft - 0.00678 Ca/sec,®C cm 

" rroxif^lno " o* 000^78 

Card No Fortran Quantity 

5 6x,dy7dz 

DX,DY,DZ are the mlnlmxim grid sizes chosen and they are 
DX - 5000.0 
DY - 5000.0 

■ "kF-T * “5^ ■ 4 ■ 

li?hen KN is the number of grids in the vertical direction 

Card No Fortran Quantity 

7 RH,kV.fi9:B'V 

These are reference horizontal and vertical eddy viscosities 

and dlffusivitles . These are calculated in section 6.3.3. 

KH » 40000.0 

KV = 10.0 

BH = 40000.0 

BV » 10.0 

Card No Fortran Quantity 

5 OmTT5'£'IT7'&l£3 T'SfifiV . EeEX , £'££Y 


45 


DELX and DELY ar« minimum grid sizes chosen « 5000 cm 
In order to determine DEEX,OEEY,EEEX,EEEY, another main 
program "COHST" has to be run twice. The Input cards needed 
for the program "CONST" are 
First time XB,A,DX,AN 
Second time YB,B.DY.BM 
Where XB « 'X boxmdary * 238000.0 cm 
A ■ a constant ■ 38,000.0 
DX - DELTA X - 5000.0 cm 
AN " Number of grids In x direction 20 
Similarly YB • 200000.0 cm 
B - 100000,0 
DY - 5000.0 cm 
BN - 20 

The output of this main program will give the constants 
C1,D and Cl.D 

First Cl,D are equal to EEEX and DEEX 
Second C1,D are equal to EEEY and DEEEY. 

For the sample problem they are equal to 
DEEX - 29333.03 
DEEY - 47503.21 
EEEX - 22947.29 
EEEY - 20957.83 

Card No Fortran Quantity 

— 5 — mThlh — 

For the sample problem the effect of wind is neglected and 
they are equal to 0.0, 0.0. If the effects wind are to be 
considered, they can be easily calculated as explained in 
Appendix A. 
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Card No 

TO— 

TAIR, Che air cemp 

Card No 

TT" 


Foreran 


n Quancicy 
TAir ^ 


30.0®C 
Fortran 


^^uantitv 


The value of DT used In the sample problem is 5 sec. In 
general It is always advised to start with a small value 
and increase it to the point when the model would not go 
unstable. 


Note: The sample problem presented is a simplified version 

of the Hutchinson Island discharge problem. 
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6,3.5, Program Ex<cutloa Procedure CNASUM III) 

This sactlon describes the procedure by which the User 
executes the NASUM III program, 

a) Input Parameter ; 

The User must first follow the steps outlined in 
the sample problem and become quite familiar with all the 
input parameters, 

b) First Run ,’ 

In order to obtain three-dimensional velocity and 
three-dimensional temperature, the main program TMAIN3 is 
executed. In TI'1AIN3 there are two tape units, One is a 
READ unit designated as Unit 7. During the first run, there 
is no need for unit 7, and unit 8 has to be provided to store 
results on a magnetic tape. 

c) Run Continuation ; 

For extending the results, the run has to be con- 
tinued, The magnetic tape which was "unit 8" in the first 
run will now be read "unit 7", for reading the previously 
stored results, Another magnetic tape is now to be provided 
as "unit 8" for storing the extended run results. The above 
procedure can be repeated until the results are obtained for 
the desired time. It is to be noted that for the first run 
IRUN - 0, and for extended runs IRUN « 1. 
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6,3,6, 

Samp la Inpuc 


Card No 

Symbol 

Value & Format 

1 

IRUN 

b b b b 0 

2 

LN 

b b 700 

3 

KSTORE 

b b b b 0 

4 

C1,CC,CP,CH,CV 

1.0, 1.0, 1.0, 1.0, 1.0 

5 

GR,FF,RR,HK 

980.0,0.00006,1.0,0.000678 

6 

DX.DY.DZ 

5000.0,5000.0,0.25 

7 

KH,KV,BH,BV 

40000.0,10.0,40000.0,10.0 

8 

DEU.DELY.DEEX, 

5000.0,5000.0,29333.03, 


DEEY,EEEX,EEEY 

47503.21,22947.29,20957.83 

9 

TAUX.TAUY 

0.0, 0.0 

10 

TAIR 

30.0 

11 

DT 

5.0 
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6.3.7. Samplt Output 

Somt of Ch« output from th« model sample run are listed 
as follows: 

a) Parameters 

b) Surface (K«I) horizontal velocities, u and v. 

c) Vertical velocity, ^2, at K»2 

d) Surface temperatures (J«l at left-hand side) 

(J-JN at right-hand side) 

I"1 at top 
1"IN at bottom 

The surface isotherms obtained after 1 hour of simula- 


tion are shown in Fig. (61 ) . 
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7.1 MAIN PROGRAM 


FOR NASUM II 


55 


1 

2 
3 
I 

5 

6 

7 

8 
9 

13 

11 

12 

13 

14 

15 

16 
17 
13 

19 

20 
21 
22 

23 

24 

25 
2S 

27 

28 

29 

30 

31 

32 

33 

34 

35 
35 
37 
33 

39 

40 

41 


C FMAIN 

C99 999 0 ****### 

C THIS PROGRAM CALCULATES THC SURFACE ELEVATIONS, CIRCULATION AND 

C TEMPERATURE OISTRIOUTION FOR THE THREE-DIMENSIONAL FREE SURFACE 

C FAR-FIELO MODEL 

t 0 0 ^ *i$0 t ^ 

PARAMETER INr ?4 , JNrl 1 ,KN =S 
REAL KH»KV 

DIMENSION U(IN, JN,KNI,V(INf JNfKNI ,U (lNviN,KN)f 
CHT(IN,JN I ,HI (IN ,JN) fETAClNf JN) tHTECINvJNI • 

CH(IN,JN,KN),C (IN, JN ,KN ) , 0 ( IN , JN , KN) ,£ ( IN , JN , KNI , 

CHXCIN,JN I ,HY( IN ,JN) ,MAR (IN, JN ) ,HTD( IN,JNI , 

CRO(IN,JN,KN) ,P(IN,JN ,KN I , HOUM ( IN , JN ) 

C,T(IN,JN,KN), TNKIN ,JN,KNI,W2(IN, JNfKNI 
C, I KIN), OTA (131 ,OHS(l3),OTAUX(13 ),DTAUY(13) 

READ 1,IRUN 
READ 1,LN 
READ l,LNl 
READ 1,M 
READ 1, INLET 
READ ItIBAY 
READ IrlHlTE 

I FORMAT (IE) 

READ ll, U,I2,I3,I<<,I5,r6 
READ 11, Jl,J2,J3,J4,J5fJ6 

II FORMAKeiCl 

READ 2, VU,TPH,TPHI,A1,A2 

READ 2,HTMIN,THT 

HEAD 2,TFLAT 

READ 2, CI,CC,Cr,CH,CV 

READ 2, GR,FF,RP 

READ 2, DX,DY,DZ 

READ 2, KH,KV 

READ 2,BH,BV 

READ 2,TI 

READ 2,(DTAUX(I ),i::l,13) 

READ 2, (OTAUY (I )vl = l,13) 

READ 2,(CTA(1 ),I=1,13) 

READ 2, (DHSa ),I = 1,13) 

DO ICO 1=1,13 
OHS(I)=OHS(I) /ICO. 


42 

43 

44 

45 

46 

47 

48 

49 
5D 
51 
5? 
53 


ICC CONTINUE 
READ 2,DT 
2 rORMAT( J 

IF(IRUN.GT.O) GO TO 4 

CALL READ2nN,JN,HAP,HI,HX,HY,0X,DY ) 

CALL iNITIAdN, JN,KN,ETA ,HT ,HI ,HTO , HTC , 1 1 ,U , V ,W , RO ,P, GR , RR , D2 ,D , 
CE,H,G,TI,T,TN1) 

TTOT=O.Q 
TT0T1=G.0 
TAUX=C.C 
TAUV=O.Q 
GO 1C 6 


54 

55 

56 


4 CCNTiNtlE 

CALL REAC'MIN ,JN,KN ,b,V ,1.' ,H I tHT ,HTD .HY ,HY ,M A H , I T A ,P , R 0 , C I , 
CCC,CM,CV ,C.i'>,DX, DY ,C2,ld,TAUy,TAUY ,1 TOT,fl,G,HTL, 1 t TTOT 1,W2 J 




57 

58 

59 

60 
61 
62 
63 
6 ^ 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 
63 

6 4 

85 

86 

8 7 
83 

89 

90 

91 

92 

93 

94 

95 

96 

97 
93 

9 9 
CO 
01 
02 

03 

04 

05 

06 

07 
03 
D9 

10 
1 1 
12 
1 2 


6 continue 

00 5 L=ltLN 

IF CTTOT.GE.THTI GO TO 1000 
TTOTsTTOT^OT 

ttsttot 

IpCM.EQ.n GO To 33 

call TIDEUN.JNfTT tLfOTfHTCfHIfHTEtlNLET^ll, 12, Jl *j 3* J4,X5tTPHl 
CtAl,A2) 

GO TO 34 

33 CALL VCL (IN «JN ,HN,U«V«H»G«O.E «TT, VO, T PH, 1 1, 12.15# Jl .J3«J4« INLET) 

34 CONTINUE 

CALL DATA (TAUX,TAUYfTT,TA ,HS ,OTAUX,OT AUYtOTA.DHS.TFLAT) 
ITTOTrTTOT 

IF (L.GT.l.OR.lTTOT.GT.O) GC TO 50 

CALL HEIGHT ( I ,IN,JN,KN ,MAP,U,V,HT ,HTD,0?*0T , OX , OY t HOU M • M , II , 
cl 2 f 13 % 14 f IS 9 16 9«il 99^3 9«i4 9w^5 9*i6) 

CALL UVVEtnN9JN9KN9U9V9H9G9DX.0Y ,DZ9W,TAUX 9TA1IY,0T ,HT , 

CHX 9 HY 9 ET A 9 P 9 HAR 9 CI 9 CC 9 CP»CH 9 CV.KH 9 KV 9 GR 9 'RR 9 Fr 9 HTD 9 R 0 ,T 9 II 9 I 29 I 3 9 
C I (4 9 I 5 9 1 6 9 1 9 2 9 J 3 9 J 4 9 J S 9 J 6 9 H ) 

CALL UVEL (IN9UN9KN.H 96 9W9HTt!l ,nx ,0Y 902 9M AR 9M 9l I9 12 9 J39 J4I 
CALL PRESI IN,JN,KN9HT09R0 9GR9P9D2 I 
CALL ETT(IN9UN9HTD9Hl9MAR9ETA) 

GO TO 80 
50 CONTINUE 

CALL HEILN(IN9JN9KN,MAR9H9G9HTD9HT,HTe, 0Z90T90X ,DY,HDUM9H9ll9l29 
CI39l4^I5 9I6 9JI 9J29U39J49JS 9U6) 

CALL UVVELNnN9jN9KN,U,V ,H ,G 9D ,E 9DX 9O Y ,0Z 9W ,T AUX 9 T AUY ,0T , 
ChT9KT09HTE9HX9HY9ETA9P9NAR9KH,KV9GR9RR9FF,CP,CC,CI,CH,CV9RQ9T9H 
C 9I29I39I49I59I69J19J29J39J49J59J69H) 

CALL WVEL«IN.JN9KN,D,E9M9HTE,nX90Y90Z9MAR9M9ll9l29J3,J4) 

CALL PRES<IN,JN,KN,HTE9R09GR,P,C2) 

CALL FTT(IN,JM^HTE,HI,MAR9ETA) 

CALL OLOUV(IN.JN9KN,U.V9H9GfO«EI 
CALL OLOHTlIfi9«N9HTE9HT09HT) 

80 CONTINUE 
5 CONTINUE 

1000 CONTINUE 

DO 15 LLrl.LNl 
TT=TT0T1*TT0T 
IF (TT.LT.THT) GO TO 3000 
IFCIBAY.EO.1) GO To 5002 
00 211 I=l9lN 

DO 211 J:i9viN 

IF{ETA<l9J).LT,0.D> GO TO 65 
GO TO 211 
65 CONTINUE 

IIII:HI ( I 9 J ) ' 

IHTHIN:hTMIN 

IF (IHI .EO .IHTMIN J GO TO 40C0 
60 TO 4n01 
4000 CONTINUE 

HTM:HTMIN»3C.48 
HI (1 9J » =HTM 

HT«l9J) :HI ( I ,J J + CTA (I ,J) 

HTC(I,J):hT(I,uJ 
HTF. n , J ) :^HT ( I 
4nol CONTINUE 





114 

115 
115 
117 
119 
119 
1Z9 
121 
122 

123 

124 

125 
12«> 
127 
123 

129 

130 

131 

132 
153 

13 4 
135 
13S 
137 
139 
139 

14 0 

141 

142 

143 


A:-30.4R 

HTMrHTNJN«30*4a 

XFICTAf It«i) .LE.A.ANO.HX(X»J).LE.HTH) CO TO 311 
CO TO 211 
311 CONTINUE 

HmtJ)sHTH«30*<*8 
HT<I,J)rETAIX(J >«HXII»J) 

HTO(ItU)-HTfI»J) 

HTEtXtU>:HTCltU) 

211 CONTINUE 
5002 CONTINUE 

IF(H.EQ*1) CO TO 330 

CALL TIOEIlNt JNtTT»LLfDT(HTOfMltHTE » INLET tIl»X2r JI,J3» J4, IS tTPHl 
C»AI,A21 
CO TO 340 

330 CALL VEL(lN,JN»KNtU(V,HrGtOyE»TT,VO»TPH*Il« I2«I5 »JlyJ3*J4f INLET) 
340 CONTINUE 

CALL OATAfTAUXy TAUT » TT , T A yH S »DT AUX y OT AUT t OT A , OHS y TFL4 T) 

XTT=TT ^ 

XTHT=THT 

IF(ITT.EQ.ITHT) GO TO 20C0 
TT0T1=TT0T1*0T 
500 CONTINUE 

CALL HEILK'(INy.JNyKN,HAR,H,GyHTOyHTyHTC»OZyDTyOX,OY,H3UMyHyIly 
C12yI3yI4yISyl6y JlyJ2,J3yJ4, J5.J6) 

GO TO 2001 
2000 CONTINUE 

CALL TIDALIINyJNyKNyCTAyHIyHTyHTDyHTEylTyT, IHITE) 

TT0T1=TT0T1*0T 
IFIIHITE.EO.I > GO TO 500 


144 

145 

146 

147 
149 

149 

150 

151 

152 

153 

154 

155 

156 

157 
159 
159 

163 
161 
162 
16 3 

164 

165 

166 
167 
169 


2001 CONTINUE 

CALL UVV£LNIINyJN,KN,U,V,HyGyOyEtOX yOYiDZyW yTAUX yTAUYyOTy 
CHTyHTDrHTEyHXyH YyETAyPyM ARyKHyKV y 6R *RR y FF y CP y CC y C I yCH « C V * RO y T y X 1 y 
C12yI3yl4 yISyI6y JlyJ2yJ3y J4y J5yJ6yM) 

CALL MVELdNyJNyKNyO yE yU y HT E yOX y OY y 02 yH AR yH y 1 1 y I 2 y J3y J4 ) 

CALL PRE SCINy JNyKNyHTEyRO yGR yPyOZ ) 

CALL ETTtlNyJNyHTEyHlyHARyETA) 

CALL OLOUVtlNyJNyKNyUyVyHyGyOyE) 

CALL OLOHTilNyJNyHTEyHTOyHT ) 

CALL TEHPlINyJUyKNyTyHyGyWyCXyOYyOZ yOTyHT ySHyBVy HARyTNlyHTOyORADy 
CRRyHSyTA yRO) 

CALL OLOT(lNyJNyKNyTyTNl) 

IS CONTINUE 
3000 CONTINUE 

CALL PRPARAICIyCHyCVyCPyCCyOXyOYyOZyOTyTAUX yTAUr yTTOTySRyFFyRRy 
CKHyKVyBHyGVyQRADyTI yTTOTl yT AyHS ) 

CALL PRETAdyJylNyJNyETA) 

CALL PRUVIIyJyK yINyJNyKN yMyC) 

CALL WW( IKyUNyKr.'yHTyHTCyETA y H y G y W yV ZyM AR y OX y OY y OZ yO Ty HX y H Y ) 

CALL PRU(INyJNyKNyMZ) 

CALL PRTEPPlIyJyK yINyUNyKNyT) 

CALL STOPEdNyJNyKNyUyVybyHIyHTyHTO yHX yHY yM AR yETA yPyROyCIy 
CCCyCHyCVyCPyDXy G V y D Z y D T y T AU X y T A U Y y T T 0 T y H y G y H T E y T y T T OT 1 y W Z ) 

STOP 

ENO 
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SUBROUTINE PROGRAMS FOR NASUM II 
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7.1.1 READ 2 


This subroutine is used by specifying 1RUN"0. First, the 
two-dimensional MAR matrix, MAR (1,J) is read in from datn 
cards (a sample problem will illustrate this in Section 6.1. 

The MAR numbering system is used for distinguishing 
between spatial differencing of the terms of the system of govern- 
ing equations in the interior of the domain of solution, on the 
boundaries, and outside the domain. Points outside the domain 
are assigned a value MAR»0, and calculations are not performed. 

The MAR matrix, as will be shown in a sample problem, is constructed 
by the user by first establishing a grid system which closely 
follows the geometry of the application site. Then the MAR 
numbering system is specified as follows : 

MAR (I,J) *0 points outside domain 
MAR (I,J) “1 upper horizontal boundary 
MAR (I,J) “2 lower horizontal boundary 
MAR (I,J) =3 left vertical boundary 
MAR (I,J) “4 right vertical boundary 

MAR (I,J) “5 through MAR (I,J) = 10 are boundary corners and 
are specified below: 

MAR (I,J) =11 interior of domain. 


MiUj. 








MAR 


5 


6 


/ > > // 

7 8 


10 
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N«xc, this subroutine reads in from data cards the two- 
dimensional matrix, HICI.J). which specifies the depth (in 
feet) below the mean water level at each grid point. Then the 
depths are converted into centimeters for calculation of the 
bottom gradients, HX(I,J) and HY(I,J) in the x and y directions, 
respectively. Note, that calculation of the bottom gradients is 
performed by using central differencing in the interior and the 
point single-sided differencing on the boundary. Again, for 
MAR“6 and MAR*8 central differencing is used, since these part- 
icular corners may be treated as being interior points. 


65 


1 

2 

3 

H 

s 

b 

7 

S 

9 

13 

11 

12 

13 

U 


c* 

c 

c 

a 


10 

99 

1 


THIS SU0R0UTIN£ READS IN DATA FOR THC APPLICATION SITE DOMAIN 
AND ACTUAL BOTTCN TOPOGRAPHY 

**«**»*»***«**•««**«*** ***««»«*«**«l****«*«»*«*««^»*«*« I 

SUBROUTINE RE AO 2 (IN ,JN »M AR * Ml ,HX ,HY •DX,OY ) 

DIMENSION HARIlNtJNltHl (IN, JNltHXIINfJNI »HY (INyJNl 

DO 10 1=1, IN 

READ 1,(MAR(I,JI,JS1,JNI 

DO 99 1=1, IN 

READ 1,(HI(1,J),J=1,JN> 

FORMAT! I 
DO 15 1=1, IN 
DO 15 J=1,JN 
HX(1,JI=30.48*HI(X,J) 


15 

IS 

CONTINUE 




15 


00 50 1=1, IN 




17 


DO 5U J=1,JN 




13 


XF(MAR(I,J).E0.0) 

60 

TO 

50 

19 


XF(MAR( I,J) .£0. 1) 

60 

TO 

31 

20 


IF(MAR(I ,J) .E0.2) 

GO 

TO 

32 

21 


IF(MAR(I,J).EC.3) 

GO 

TO 

33 

22 


IF(HAR(I,U),£O.MI 

60 

TO 

3** 

23 


1F(MAR(I,J).E0.5> 

GO 

TO 

35 

24 


XFiMARd, J).E0.6) 

GO 

TO 

36 

25 


IF(HAR(I,J).E0,7) 

GO 

TO 

37 

25 


1F(HAR(I,U) .EQ.6) 

GO 

TO 

38 

27 


1F(MAR( I ,J) .EC. 9) 

GO 

TO 

39 


23 ir(HAR(I,J).EO* 10) 60 TO MO 

29 HX(I,J)= (HI(I*l,J)-HI(I-l ,J»)/ (?.*DX1 

50 HY(I,.J)=(HI(X,J«n-HI(I,J-l ) )/(2,Y0Y) 

31 CO TO 50 

32 31 CONTINUE 

33 HX(I,U)=(HI (I«l ,J)>HI(I*1 ,J ) 1/ (2.*DX) 

34 HY(I,J) = (3*HI(I ,J)«HI(I, J*2 )-M*Hl(I,J>ll )/(2«DY) 

35 60 TO 50 

35 32 CONTINUE 

37 HX(I,J)=(HI(I41 ,J)>HI(I-l,Jn/(2.*0X > 

33 HY(I,J)=(M«HI (I ,J4l)-3«HI (I ,J)-H1(I ,J«2 l)/(2«0Y) 

39 60 TO 50 

43 33 CONTINUE 

4 1 HX(I,J)=(4*HI(I«l,J)>3>^HI(I,J)*HI(I<»2f J) )/(2«0X) 

4 2 HY(I,J) = (HT(I,J«1)-HI(I,J>1 n/(2.«DY) 

43 60 TO 50 

44 34 CONTINUE 

45 HX(I,JI=( 3*HI (I ,J)«HX(I-2,J )-M*HI(I -l,J) »/(24DX ) 

4 5 HY(I,J)=(HI(I,J4l)>HI(I.J-I ))/ (2.<*0Y) 

47 60 TO 50 

48 35 CONTINUE 

49 HXII,J)=(M*HI(I*l,J)-3*HI(I,J>-HI(l42,J))/(2*DX) 

5 0 HY(I,J)=(3*HH1 ,J |4HI(I , J-2 )-M *HI ( 1 , J-U ) / ( 2*DY ) 

51 GO TO 50 

52 36 CONTINUE 

53 HX(1,J)=(HI(I41 ,J)-HIII-l ,J ))/ 12 .♦OX) 

54 HY(I,JI=(HI(I,J*1)-HI(I,J-I))/(2.’«’DY) 

55 GO TO bO 

55 37 CONTINUE 
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57 


HX(X|Jt:f 4*HX (X «ltJl-3*HlCI tJI'HX (X«2»JI l/(2*0X> 

55 


HYtXiJIsUwHI IX »J«n*3*HXIX (JI'HIII »U«2 1 1 / 1 2«DV I 

59 


60 TO SO 

59 

38 

CONTINUE 

51 


tlXIXf J)S tHIIX«l*J}>HIfI*l(jn/C2.*0X) 

52 


HY(X,J)=IM1(X |J«U>HX(XtJ-l)l/<2*<*0T) 

53 


60 TO SO 

55 

39 

CONTINUE 

55 


NXIXtJ)St3*H! IX *J }«HX IX«2 , J l-4*HX 11 -ItJI 1/12*0X1 

55 


HYIX»Ji:i54>HX II tU«ll-3»HllX tUl-HIIX »J«2 1 1 / 1 2*0Y 1 

57 


60 TO SO 

55 

50 

CONTINUE 

59 


HXI X » J1 = 1 3*HX IX tJl*HXIX>2 »Jl-4*HIIX-i»Jl 1/1 2*0X1 

79 


HYIX,JUI3*HX IX tJl«HllltJ-2 l•4*HXIX ,J-1 1 1 / 1 2*0Y 1 

71 

SO 

CONTINUE 

72 


RETURN 

73 


ENO 










'WW 
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7.1.2 INITIA 

This subroutine is used by specifying IRUN "0. The initial 
values of n, H, u, v, n, T, P and p are specified as was out- 
lined in section 4.1. 
n - 0 
u " 0 
V ■ 0 

■ 0 

*T « - const. 

P " Pq + PRHc 

“ ^unif 

In terms of the symbols used in the program, we have: 

ETA(I.J) - 0 
HT(I.J) - HI(I.J) 

HTD(I,J) - HI(I.J) 

HTE(I.J) - HI(I,J) 

U(I,J.K) - 0 
H(I.J.K) - 0 
D(I,J.K) - 0 
V(I.J.K) - 0 
G(I.J.K) - 0 
E(I,J,K) - 0 
W(I,J,K) - 0 
T(I,J,K) »TI 
TN1(I,J,K)-TI 


RO(I,J,K) »RR 

P ( I . J , K) »RR*GR*HT ( I , J ) '4C- 1 ) *DZ 

*NOTE: The energy equation is only coupled to the system of governing 

equations after the IR data base is inputted at 10:10pm for 4/15/75 
into the program, and therefore an initial value of T is not actually 
required. However, if no IR data is available an isothermal bay may 
be assumed initially and the coupling to the energy equation initially. 
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t 

2 

S 

1 

s 

» 

T 

8 

9 


C 

c 

c< 

c 

c« 


THIS PROCrAM iNITIALIZrS VARIABLES TOR CONSTANT OCNSITV NOOCL 


THIS SUPRCUTPC IMTIALIZCS THE VARIABLES FOR THE VARIABLE DENSITY 

SUBROUTINE INITI A UN , JN ,KN, ET A ,HT »H1 ,HTD ,HTE » II * U»V»M »RO, P»CR i RR • 
COZtOte*H*C»TItT»TNn 

OINENSION ETAII r;,JN),HT( IN| JM,HI(IN,JNI ,HTD(IN*JN)*NTEI IN* Jt. >« 
CimNI,U(lN*JN,KN »*V(IN*UN,KN) ,D(1N *JN*KN l*E(IN»JN*KNI*U(IN,JN,KN) 


13 


C *R0( INtUNfKNI »H IIN*JN|KN ) *G( IN , 

II 


CTNlCINt JN,KNI 

12 


DO 2C isitlN 

13 


DO 20 JSlfJN 

t<l 


ETAII»J)70.0 

IS 


HT(I*J):HI(I* J| «ETA(I,J| 

IS 


HTO(ItJt:HT(I»J) 

17 


MTEtl»J|:HT(ltJ> 

19 

20 

CONTINUE 

19 


DO e IsliIN 

20 


DO 6 JsifJN 

21 


DO 8 K:1,KN 

22 


UIItJ*K|:C* 

23 


MfI*J*K):C* 

29 


D(I»J(K):0. 

2S 


V(Zf JfKiro. 

26 


G 1 1 f «if K ) *0« 

27 


Cll*J»KI=0. 

29 


U(I,J*K)rO. 

29 


T(I*J,N)=T1 

JO 


TNK I*J*KUTI 

JI 


RO(I,J»K|rl.029931-.C0Q020*T (I, 

J2 

8 

CONTINUE 

33 


DO 2S ISI.IN 

39 


DO 25 jri*JN 

35 


DO 21 K=1*KN 

36 


RR=RO(If J*K) 

37 


PIItJtK}=CR«HT(I»J}*RR«(K-l )*02 

33 

21 

CONTINUE 

39 

25 

CONTINUE 

90 


RETURN 

91 


ENU 
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7.1.3 READ I ] 

i 

This subroutine Is used by specifying IRUN"1. The results I 

j 

of a previous run are read In from a magnetic tape for the 
purpose of running the program over a long period of time In 
segments. The system variables n, H. u, v, T. P and o 
are read In as well as physical and numerical parameters, as 
follows : 

" depths at each grid points (in cm) 

HX(I,J) “ bottom gradient in x-direction 

HY(I,J) ■ bottom gradient in y-direction 

MAR(I,J) - domain numbering system 

Cl, CC, CH, CV, CP ■ constants always » 1 

DX » horizontal grid size in z-direction (in cm) 

DY » horizontal grid size in y-direction (in cm) 

DZ ■ vertical grid size - AZ/HI(I,J) 

*DT “ time step (in seconds) 

TAUX * wind stress component in x-direction (dynes/cm*’') 

TAUY ** wind stress component in y-direction (dynes /cm") 

TTOT “ total run time without energy eouation (in seconds) 

TTOTl » total run time with energy equation (in seconds) 

* NOTE: The determination of At is presented in the sample problem 

in section 6.1. 



1 c 

3 C THIS SUBROUTINE READS XN(FROH MAGNETIC TAPE) VALUES FOR THE VARIABLE. 

R C AND PHYSICAL AND NUMERICAL PARAMETERS FROM A PREVIOUS COMPUTER RUN 

5 C FOR THE VARIABLE DENSITY MODEL 

6 C^******************^^•*****0*m*m****^*****0**0******^^*^0**^>*i,•w»^*^**i^**9*■‘ 

7 SUGROUTINE READ 1 (IN,UN,KN»U »V,W,HI»HT ,HTD(HX •HY»MAR»ETA»P»RO,CIt 

I CCC«CH,CVtCPtDX»OY*D 2 tOT»TAUX,TAUV,rTOT,H»G*HTE|T,TTOTl»UZ) 

9 DIMENSION U(I N, JN »KN ) , V ( IN* JN ,KN I IIN, JN f KN I tP ( XN» JNt KN) , 

ID CHI(XK»UNItHTIIN»JN) tHTO(INfiN),HX(INyJN> »HY (IN,JN) yMARtlNy JNl, 

11 CCTAIIN,JN)tROaNtJNyKN)»H( 2 N,UN»KNt »G(IN , JN »KN) , HTEIINy JNl 

12 C|T(TNyJN(KN)fWZ(INyJN»KN) 

13 READ t 7 > mU(I,JyK},K:itKN),j:ltJN)yZrl,XNt, 

14 CU(V(IyJ«K)tK:l»KNUJ:itUNI yIslyIN) f 

15 Ci ((WfIyJ,K),K=UKNItJriyUNIylrlyIN) „ 

IS Cl l(H(ItJfK)»K:itKNI«J=I* JN» ylslylN) y 

17 Cl I (GlIyJyKIyKSlyKN) y U= I y UN ) y I S 1 y I N ) y 

IS Cl IIPlIyUyKtyKSlyKMIyUSlyJM) yISlyIN) y 

19 Cl I IROlIy JyKIyKSlyKN lyJSlyJN )yl = l yIN )y 

20 Cl IHTDIIy JlyUSlyUNlylSlyINIy 

21 Cl IHTEIlyUlyUSly JN)yI=lyIN), 

22 CI|HX(IyUlyU=lyUN)yl=lyIN)y 

23 Cl IHX{IyJ)yJ=IyUN)yI=lyIN i y 

24 CIIHYIIyUlyJ=lyUN)yI=lyIN>y 

25 Cl IHARII, JlyJ=lyUN)yI=l ,IN)y 

2 S CIIHT(lyUlyJ=lyJN)yI=lyIN)y 

27 Cl lETAlIy JlyU=lyJNIyI=lyIN)y 

28 Cl (IT(lyUyK)yK=l,KN)yJ=lyJNI yI=lyIN) y 

29 Cl I |WZIIyUyKlyM=lyKN>yJ=lyJN >yI = lyIN )y 

30 CCIyCCyCHyCVyCPyCXyOY.DZyDTyTAUXyT AUYyTTOT iTTOTl 

31 RETURN 

32 END 
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7.1.4 TIDE (M-2) 

This subroutine specifies the tide height at the ocean-bay 
opening as a function of time (i.e. TTOT or TTOT + TTOTl) in ttie 
form: 

(t) - + A 2 cos fi)(t+t(j)) 

or in terms of the program symbols: 

ETA - A1 + A2 *cos ( (TT + TPH1)*(6. 23/12. 15 *3600)) 
where A1 and A2 are read in from data cards , TT ■ TTOT or 
•TTOT + TTOTl and TPHl is computed as will be illustrated in the 
sample problem. 

This is the actual tidal condition case for tidal flux at 
the ocean-bay opening. Note, this subroutine is not used when 
the current velocity is specified at this open boundary (i.e. 

M“i for specifying V^(t)). It is further pointed out that 
the user must start the computer run at t»0 (TTOT»0) with a flat 
surface). ETA(I,J)=0, and ETA^O at the open boundary. Otherwise, 
the large step in surface height at the open boundary will 
result in numerical instability, since the governing equations do 
not adjust sufficiently fast to yield a compatible and realistic 
situation between surface heights and currents. Since, the dom.ain 
of solution is assumed to be "still" (i.e. u*v=fi*o) initially 
for the case of not having an adequate initial data base for spec- 
ifying currents, the above noted specification of ETA(I,J) at 
t=0 must be followed to insure numerical stability. 



1 c«****<»*»********'*>» **•*****'»*«'*«**'*•«'* to* *•*«*•«•**•*«* ******•***<*«•***««*»< 

2 C THIS SUnROUTiKE CALCULATES THE SURFACE ELEVATION AT THE OCEAN-BAY 

5 C OPENING AT EVERY TIME STEP 

4 C *•***<*•*<**•>***»»«>»**••***•******>»•«*********«****«*** **»**«*«>***«**»***«x 

5 SUBROUTINE TIDE I IN » JN , T TO T , L ,DT « HT D |HI «HT C t INLET ,I1»12,J1 ,J3,J4 

6 CtlS.TPHl »A1,A2) 

7 OINENSION HTDIIN,JN),HI(IN, JN),HTE(IN,JNI 

S ETA=Al«A2«COSIiTTOT«TPHl ) * i 6 . 26/ ( 12 .15 *36 00 .0 1 1) 

9 XFdNLCT.CT.l ) CO TO 2000 

10 00 40 l=IltI2 

11 1TT0T=TT0T 

12 IF<L.CT. l.OR.ITTOT.CT.G) GO TO 5 

13 HTO( I,J1)=ETA«HI(I,J1) 

14 GO TO 40 

15 5 HT£( I»jn=ETA«HI(I»Jl) 

16 40 CONTINUE 

17 GO TO 1000 

18 2C00 CONTINUE 

19 DO 50 J=J3»U4 

20 ITTOT=TTOT 

21 inL.GT. l.OR.ITTOT.GT.O) GO TO IS 

22 MT0(15»J)=ETA«HI(I5.J) 

23 GO TO 50 

24 15 HTE(I5»J)rETA«HI(IS»U) 

25 SO CONTINUE 

25 ICOO CONTINUE 

27 RETURN 

23 END 
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7.1.5 VEL(M»1) 

This subroutine specifies the current velocity at the ocean- 
bay opening as a function of time ( ie TTOT or TTOT + TTOTl) 
in the form: 

Vq (t) • cos 0 ) (t + ti<>) or; 

Uq (t) ■ cos 0 ) (t + t(j>) depending on which axis the open 

boundary is aligned. 

Then, in terms of the program symbols: 

Vo - Vo*cos ((TT + TPH) * (6 . 24) /(12 . 15 * 3600)) 

Where Vo is computed as an estimate and TPH is computed as will 
be illustrated in the sample problem. TT = TTOT or TT = TTOT + 
TTOTl. Similar!,', Vo (t) » 0 at t =* 0 is recommended along with 
an initially 'flat' surface to insure compatibility between 
the surface heights and the velocity field. 


o n o o 


THIS SUBROUTINE CALCULATES THE CURRENT VELOCITY AT THE OCEAN-BAY 
OPINING AT EVERY TIPE STEP 


SUbROUTINE VEU INtJNfKNtUfV tHvCfOfE «TT«VO vtPHfll tlZvISf JltJSf J4t 
C INLET) 

DIMENSION U(XNt JNfKN) ,V( IN» JNvKN » fH flN t JN »KN ) vGf IN t JN «KN ) 
CfOdNfJNfKN) f E ( IN tUN |KNI 
KNUHN-l 
DO 1000 KrliKNl 
XF(INLET.GT.l) CO TO 2000 
DO 100 I^Ilfl? 

VC It Jit K )rVO«COS( CTT4TPH)«(( *26/(12 *IS «3 6 00 *0 I ) I 
CiCit JlfK )=V(I tJ IfK) 

E(IttUfK)=G(I tJlfH) 

ICC CONTINUE 

GO TO 1000 
2000 CONTINUE 

DO 200 j::J3tJ4 

U( 15 t JfK )=VO*COSC (TT«TPM ) *( 6 • 26/ ( 12 .15 b 00 *0 ) ) ) 

HCIStJfK)=U(I5t JfK) 

OCI5tJtK)=:HCI5t J,K) 

2CC CONTINUE 
1000 CONTINUE 
RETURN 
END 
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7.1.6 TIDAL 

This subroutine specifies T(I,J,K) as an initial temperature 
distribution which is constructed by the user from an IR-Data 
Base. This will be illustrated in the sample problem. Spec- 
ification of an initial surface as constructed from an adequate 
tide data base may be inputted into the model but is left 
optional for the user. 
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13 
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32 


C****4i«***9*«***«i9«*«i**99*»*99*»**«9**»»«*»*9 9«»*»»*#»***M*9'»**«4i*9«**«>»«. 

C THIS SUBROUTINc READS IN THE INITIAL VALUES FOR THE TEHPERATURC 

C DISTRIBUTION AND THE INITIAL VALUES FOR THE SURFACE 

C ELEVATIONS(OPTICNAL) 

C9*0*«*»*»**« ***** *****«*»*****9 *•*•*»** t********************************' 

SUBROUTINE TIDAL IIN »JN «KN ,C T A ,HI , MT *HT D » HTE * 1 1 » T « IHIT E I 
DIMENSION IKINTtETAdN, JNt tHI (IN«JN),HT IIN»JN)»HTD(IN, JNt» 
CHTL(IN»JN),TIIN,JN»KNI 
DO 100 I=1»IN 
READ S»mi«J,n,Jsl»JN) 

5 FORMAT! I 

ICO CONTINUE 

DO 101 1=1(IN 
DO 101 JrifJN 
DO 102 K=2»KN 
TCIf JtK)=T(I,Jt 1) 

102 CONTINUE 

101 CONTINUE 

IFdHITE .E0*1I GO TO 25 

DO 10 1=1»IN 

READ l*Il(I)t(ETA(I,J),J = l,JN) 

1 F0RHAT(I3,11F7.2) 

1C CONTINUE 

DO 20 1=1, IN 
DO 20 J=1,JN 

HT(IyJ)=ETA(I ,J I«H1 II,J) 

HTO( I,JI =HTII ,J ) 

HTE(I,JI=HT(I,J t 
20 CONTINUE 

25 CONTINUE 

RETURN 
END 
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7.1.7 DATA 

This subroutine specifies Che wind stresses, equilibrium 
temperature and the surface heat transfer coefficient every 
hour. However, these physical parameters must be calculated by 
Che user as will be shown in the sample problem ^ Note, that this 
subroutine has been programmed for a maximum computer run of 12 
hours. However, this subroutine can be used for less chan 12 hours, 
or it may be easily modified by Che user for computer runs in 
excess of 12 hours, by merely reading in values of these physical 
parameters for the desired increase in number of hours , 
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1 

2 C THIS SUBROUTXKC READS IN HOURLY VALUES OF SURFACE WIND SHE^^ STRESS* 

3 C EQUILIBRIUM TEMFERATURE* AND SURFACE HEAT TRANSFER COEFFICIENT 

5 SUBROUTINE D A TA I T AUX ,TAUY »TTOT »TA,HS lOTAUXtOTAUV * DT A^ DNS * TFL A Tl 

6 DIMENSION OTA UX I m tDT AU Y <1 3 1 * OT A CX 3 ) |0HS (1 3 ) 

7 TTOTrlTOT-TFLAT 

a IF(TT0T«LT*0.0) CO TO SO 

9 IFITTOT.GE.O.r. ANO.TTOT. lt. 3600. 01 GO TO I 

XO IF(TTOT.LT.72CO.ANO.TTOT.GE.3600«.> GO TO 2 

11 IFCTT0T.LT.1080C.AND.TT0T.GE.7200. > GO TO 3 

12 IFdTOT.LT.H^Or.AND.TTOT.Gc.IOBOO. I 60 TO R 

13 IFCTT0T.lt. 18C0D. AND. TTOT. 6E. HRUO. ) GO TO 5 

n IF(TTOT.LT.21C3C.AND.TTOT.6E.16COO. 1 GO TO 6 

15 IFCTTOT.LT.2520C.AND.TT0T.GE.216CC. 1 CO TO 7 

16 IFCTTOT. LT. 20800. AND. TTOT. GE. 252C0. > GO TO 8 

17 IFCTTOT. LT.32«*00. AND. TTOT.GF. 28800. ) 60 TO 9 

18 tfcTTOT. LT.36C0C. AND. TTOT.OE. 32400. ) GO TO 10 

19 IFnT0T.LT.3963C.AND.TT0T.GE.:6000. I GO TO 11 

20 IF(TTOT.LT.43?OC.AND.TTOT.CE.396CO. I GO TO 12 

21 XF(TTOT.LT.460CO.AND.TTOT.GE.432CC. ) GO TO 13 

22 1 CONTINUE 

23 TAUXrOTAUXCl J 

24 TAUYrOTAUYCl) 

25 TAs^DTACn 

26 HSrOHSIl ) 

27 60 TO 50 

23 2 CONTINUE 

29 TAUXrDTAUX(2l 

3D TAUY=DTAUY(2) 

31 TArOTAI2) 

32 HS=DHS(2) 

33 60 TO SO 

34 3 CONTINUE 

35 TAUXrDTAUX(3> 

36 TAUY=DTAUY<3) 

37 TArDTA(3l 

3S HS=DHS(3 I 

39 GO TO 50 

40 4 CONTINUE 

41 TAUX=DTAUX<4 I 

42 TAUY=0TAUY14 ) 

43 TA=0TAC4 I 

44 HS=DHSC4 ) 

45 GO TO 50 

46 5 CONTINUE 

47 TAUXrDTAUX(S) 

43 TAUY=OTA UY(5 ) 

49 TA=DTAC5) 

53 HSiDHSCS I 

51 CO TO 50 

5? 6 CONTINUE 

53 TAUX-DTAUX<6 ) 

54 TAUY::0TA UYC 6 ) 

55 TA=rTA(6 ) 

56 H5-DMSC6 ) 


r 

I 
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i 


57 


GO TO 50 

55 

7 

CONTINUE 

59 


TAUXSDTAUX(7) 

63 


TAUY=OTAUYI7l 

61 


TASOTAtT) 

62 


HSSOHSIT) 

63 


60 TO 50 

65 

8 

CONTINUE 

65 


TAUX=0TAUX(8I 

66 


TAUY=DTAUY(B) 

67 


TA=DTA(8I 

6B 


HSsOHS(8 1 

69 


60 TO SO 

73 

9 

CONTINUE 

71 


TAUXroTAUXI9) 

72 


TAUY=0TAUY(9) 

73 


TAS0TAC9I 

75 


HS=0HS(9) 

75 


60 TO SO 

76 

10 

CONTINUE 

77 


TAUX=DTAUX(10 ) 

75 


TAUYSOTAUY(IO) 

79 


TAsOTAt 10) 

83 


HSSOHS(IO) 

81 


60 TO 50 

82 

11 

CONTINUE 

83 


TAUX=0TAUXI11 1 

85 


TAUY=OTAUY( 11) 

85 


TASOTACll) 

86 


HSSOHS(II) 

87 


GO TO SO 

88 

12 

CONTINUE 

89 


TAUX=0TAUX(12 ) 

93 


TAUY=0TAUY(12) 

91 


TA:0TA(12) 

92 


HS=DHS(12) 

93 


GO TO 50 

94 

13 

CONTINUE 

95 


TAUX=DTAUXI13) 

96 


TAUY=OTAUY( 13 ) 

97 


TA=0TA(13) 

93 


HS=0HS(13) 

99 

50 

CONTINUE 

133 


TT0T=TT0T4TFLAT 

131 


RETURN 

132 


END 


preceding 


elank 


filmed 


I 


Kk.'- 


vi.rrv' 


• w 
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7.1.8 HEIGHT 

This subrouclne calculates H(x,y) the depth contour with 
respect to the free surface, at t ■ At (■HTD(I,J)> by forward 
differencing the surface height equation in time with respect 
to the initial depth contour matrix H(x,y) at t»0 (-HT(I,J)). 
Mote, this subroutine is used only for the first time cycle. 

The integration in this subroutine is performed by applying 
Simpson's Rule.‘>^ The general inlet and outlet conditions are 
specified by reading in from data cards parameters which impose 
the location of the inlet, either on the upper horizontal 
boundary or on the left vertical boundary of the grid system. 
However, this subroutine can be easily modified by the user for 
having an inlet on the lower horizontal boundary or on tne right 
vertical boundary. The only change required is respecifying 
MAR(I,J) corresponding to the inlet location and reading in from 
data cards the values of I and J which properly locate the inlet. 

The derivatives in the integral are obtained by central 
differencing in space for interior points, including MAR *6 
and MAR- 8 . Three point single sided differencing is performed 
on the boundaries. These differert schemes are given in Volume I. 

jF(x)dx ” Axl 1 / 3 F(x 3 _) + 4 / 3 F(x 2 > + 2/3 F (X 3 ) 

A +4/3 F(x^) + 1/3F (X 5 )] 


* For KM » 5: 
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1 c 

S C THIS SUBRCUTXNC CALCULATES THE TOTAL DEPTH AT EACH X-V LOCATION 

B C IN THE DOMAIN FOR THE FIRST TIME STEP USING A FORWAR3 OIFFERCNCXNG 

S C SCI.CHE IN TIME 

7 SUBROUTINE HEIGHT (I »J,K tlNfUN #KN |MAR »UtV tHT f HTD« 07 tCT f OX tOV (HOUM , 

I CM(Iifr2(I3tI^tIS|l6|Jl|J2tJ3 t<AR t <J5 f Jb I 

9 OIHCNSXON HAR(INtJNI,U(INtJN,KN|,V(INtJN«KNItHT( INtJMI(HTOIIN|JNI, 

19 CHOUMlINtJNI 

11 KNHUKN-I 

12 00 SO IsltlN 

IS 00 SO JSltJN 

IH HOUMdtJisO.O 

15 00 6C KS1,KN 

16 XF(HAR(1(J).E0.C) CO TO SC 

IT XFIMARlIfJI.EO. 11) 60 TO 11 

19 XF(HAR(l,J).eO.S) GO TO 12 

19 XF(MAR(XfU).EO*S) GO TO 19 

20 XFIMARt I,J).EO.Z) CO TO 13 

21 XFIHAR(I,J).EO.(<) GO TO 14 

22 IFIHARt X tU)*EO. 1 ) GO TO 2C 

23 IF(MAR( I ,J) .EO. 7) GO TO 15 

24 XF(HAR(I«J) .E0.9) GO TO 16 

25 XFIMARIXtUI^CO.lOl GO TO 17 

26 XF(HAR(I,J).E0.6) GO TO It 

27 XF(MAR(XtJ)*C0.P) CO TO 11 

29 11 01HUX = (HTI1«1,J l»U(I-»l,U,K)-HTII>l,U)*UIX>l,J(K) )/(2.*0X) 

29 DlHVY=IHT(ItJ«l )*V(I,J*ltK)-HTII,J>l>*V<IyJ-ltK) l/t2.*DVI 

50 GO TO 24 

51 12 CONTINUE 

32 IFd.EO. I5.AN0. J*GE*J3.AN0. J.LE.J4. ANO.H.GT.l) 60 TO 50 

S3 0UIUXS(4*HT(X«1 •J)*U(I«lyJ»K)>3«HTI I» J) *U ( I , J,K) -HT( I« 2y J)* 

34 cuii«:tUyKn/( ?.*ux) 

35 01HVY=IHT(ItJ«l )9VII tJ«lyK)>KT(I,J*l)4V(ltJ*lyK) )/(2.4DY) 

36 GO TO 24 

37 14 CONTINUE 

39 01HUX=(34HTfX yJ )*U(I yJyK )«HT(I-2y J) *U(I*2yJyK|-4*HT(I>ly J)* 

39 CU(I*lyJyK))/(2t*0X) 

4 0 DlHVY=(HT(IyJ«l)*Vn yJ«lyK>-HTIIyJ-l)*tf <1 yJ'lyKI )/l2*<<0Y) 

41 GO TO 24 

42 13 CONTINUE 

4 3 01HUX=(HT(I«lyJ )*U(I«lyJyK)>HTIl'ly J)*U(I*lyJyK) )/(2.40X) 

44 0lHVYr(4*HTI I y J « 1 l*V ( I y J « 1 y K )-3*HT ( I y J) *V I I y J yK ) -HT ( I y U« 2 ) * 

45 CVIIyJ«2yKn/C.*0Y) 

46 CO TO 24 

47 20 CONTINUE 

49 IF(J.EO. Jl.ANO.I.GE.Il.ANO. I.LE.I2. AND.H.GT.I) CO TO 50 

4 9 01HUX=(HT(I«lyJ >«U(I«lyJyK)-HT(I-lyJ)*U(I-l yJyK) )/(2.*0X) 

50 01HVY=C3*MT(I yj )*V(1 yJyK )«HT (IyJ-21 «V I I y J <>2 y K > -4 »HT ( I y J- 1 ) « 

51 CVI ly J-ly Xn/( 2.40Y) 

52 GO TO 24 

53 15 CONTINUE 

54 01HUX=(|4*HT(I«lyJ t»U(I«l y Jy K )>3 4HT( IyJ)«U (I yJyK) -MT(I«2y J)« 

55 CU(l«2yJyK)t/( Z.^^DX) 

56 DlHWYrCM *HT( I y J «1 )*V«1 y J ♦ 1 y K ) -S*HT I T y J » *V « I y J yK ) -NT « I y J* 2 ) ’!■ 
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ST 


CVIX»U«2,K)>/( 2,90Y) 

Sf( 


CO TO 29 

S9 

19 

COUTINue 

S3 


01HUX = (9 4HTII«ltJ)9Ull«l,J,K)-3*HT( I»J)«U(XtJ»KI 'HTIMZi Jl* 

SI 


CU1I«2(JtKI}/(2.90X> 

ST 


0lMVYSI3 9NTII,J)«VfI (JiK )«HT(I»J>?l*VfIt«i'ZtKI>9*HT(I»J>l)9 

S3 


CV(l»J-l|Ktl/| 2.90YI 

ss 


60 TO 29 

S5 

16 

CONTINUE 

SS 


01HUX=(3*HT(1 tJ 19U(I »J|K l«HT(l-Z t J) *U ( X -2 « J t K t -9 *HT ( I -1 r J 1 * 

ST 


CUIX-XtUfNII/(2.*0XI 

S3 


DlHVV = t9*MTf X *J«1I*VII ,J«1,KI-39HTI I,J)9V (X »J|K) -HTiXt J*2>* 

S9 


CVIXt J«2»KII/I 2.*0Y) 

13 


GO TO 29 

T1 

17 

CONTINUE 

TZ 


OlHUXZISfHTlX |J >9UII iJiK » «HT ( I -2 i J) *U II -2 t K ) -9 *HT (I -1 , JM> 

T3 


CUII-l,J,KM/< 2. *0X1 

79 


0XHVYrf3*NTIX fJ l*V(I ,J,K 1 «HT 1 1 » J-2 ) *V (I t J-2 t K 1 -9*HT I X * J-1 > « 

75 


CV(I,J>1,KI)/I 2.*0V) 

7S 


GO TO 29 

77 

29 

CONTINUE 

78 

c • • 

..SXNPSON*S RULE IS USED FOR INTEGRATION 

79 


IFIK.EG. l.OR.H.fO.S) GO TO 101 

83 


XFIK.EQ.2.9R.K.C0.9) GO TO 102 

81 


HOUMIXtJUl (OlHUX«DlHVY|aOZ«l2*/3» > >«HOUHf I , J> 

82 


GO TO 103 

83 

101 

M0UMIX,J)=l(0 lHUX«01HVVI«0Z/3. MHOUMII, J) 

89 


CO TO 103 

85 

102 

H0UNII»JI=ltDlHUX«01HVY)*0Z*(9./3.) )«HOUM(I» J) 

8S 

103 

CONTINUE 

57 


60 CONTINUE 

83 


HT0IX»J)=HT(1»J l-HOUPdtJ) 

89 

SO 

CONTINUE 

93 


RETURN 

91 


END 


iNTE4mc;jW mcA 


86 


7.1.9 HEILN 

This subroutine calculates H at time level n-il (■HTE(I,J)) 
from H at time level n("HTD(I,J)) and H at time level n-l("HT(I ,J)) 
by using central differencing in time. 

Volume I gives the detailed finite difference scheme used by 
this subroutine for solving the surface height equation. The 
integration, once again, is performed by using Simpson's Rule. 

The general inlet and outlet conditions are incorporated as was 
just discussed for subroutine HEIGHT. 
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C 

C****«*****t****«*****»*9*****«***94>*i»*****»4M»*»«*«4>«4ii^4*4>»ti4i***i)i**i»*«*44<>'(' 

C THIS SUBROUTINE CALCULATES THE TOTAL DEPTH AT EACH X>Y LOCATION 

C IN THE DOMAIN P OR THE SECOND TIME STEP AND THEREAFTER USING A CENTRAL 

C DIFFERENCING SCHEME IN TIME 

C****«.«*«*******«i*«*4M»«*******»*«*i»***9*»«»***9**99*«>9**4>*A*«>994>****9*>t>*««-ik 

SUBROUTINE HEIL N (IN »JN fKN,K AR ,') » V tHT»HTO»HT E »02 (DT »DX , 3Y , HOUH» H 
C(Il|I2,13|I4t IS tI6,Jl,«l2,J3*J4,JS,J6) 

OIHSNSION HAR(lN(JNItU(IN»JN,KN),V(IN»JN»KN)tHT(INtJN}»HYD(lNtJN)( 
CHOUHIXNt JN)»HTE (ZNiJNI 
KNHlrKN*! 

DO 50 l:i(IN 
DO SO JSlfJN 
HOUM(I»J>:0.0 
DO 60 K:1,KN 

IF(HAR(T iJi .EO. 0) GO TO SO 
IF(MAR(I,J).EO. Ill CO TO 11 
IF(HAR(ItJI.E0.3l GO TO 12 
IFIHARd.JI.EO.S) GO TO 19 
IF(MAR(I ,JI.EQ.2) 60 TO 13 
IF(MAR( I tJ) *E0. 1) 60 TO 20 
IF«MAR(1»J).EC.4> 60 TO 14 
IF(MAR( I , J) .EO. 7) 60 TO IS 
IF(MAR(I(J).E0.9> 60 TO 16 
IF(KAR(I,JI.E0.1QI GO TO 17 
IFtMARd ,J) .EO.ei GO TO 11 
ir(MACdtJI*E0.8l 50 TO 11 

11 01HUX=(HT(I*1,J l*UII«l,J«K)-HTd-ltJI*Ud>ltJ«K) l/(2.*0X) 
DlHVY=(HT(I,J*l >4Vd,j4l»K}-HT(I,J>l)*Vd|J-l»K) l/(2.«DYI 
GO TO 24 

12 CONTINUE 

IFd.E0.I5.AND. J.GE.J3.AN0. J*L£*J4.AND.H.6T.1) GO TO 50 
01HUXr|4l>HT(l4l «J l«U(l« 1 « J»K l-3*HT( 1 1 J I *U (I » J fK I -HTd * 2| Jl • 

CUd«2»Jt K))/( 2. *0X1 

01HVY = (HT(I,J41 )•V(I,J♦l,KI-HTd,J-ll'»Y(I tJ-l,KI l/(2.4DY) 

GO TO 24 

14 CONTINUE 

01HUX = (34HT(I ,J M«U(I,J,K)4HTd-2, Jl*Ud-2,J,K)-4*HT(I-l, Jl* 
CUd-l»J*K))/(2.«0X) 

01HVY=(HT(I |J«1 )9V(I ,J«IyK)-HT(I t J>l)*Vd 1/(2.4>0YI 

GO TO 24 

13 CONTINUE 

01HUX = (HT(I«1»J )*U(I41 |J tK)-HTd«ltU)*Ud-l tU»K| l/d.*OX) 

D1HVY=(4 *HTd ,J n l•*>V(I,J-»l,K)-3*HTd,J}*V (I «J,K) -HTII, J*2I* 

CVd.J42, K))/( 2.4'OY) 

GO TO 24 
20 CONTINUE 

IFCJ.EQ. J1.AND.I.GE.I1.AN0.I.LE.I2.AN0.H.GT.1) 60 TO 50 
D1HUX=(HT(I4 1.J )♦U{I4l,J,KI-HTd-l, J)*Ud-l ,U,K I I / ( 2. ♦ DX I 
D1HVY=(3*HT{I ,J I^VdyJ.K J 4HT (I , J~2 J *V ( I , J -2 ,K )-*4*‘HT ( I , J- 1 ) * 

CVd, J-l, Kl) /( 2 . 'aOYI 
GO TO 24 

15 CONTINUE 

01HUX=(4 *HT(l4l ,J >*U(l4l, J, KI-JfHTdyJlOU (I , J»K) -HT(I42, Jl -A 
CUd42tUtMn/(2.«DX> 

D1HVY=(4 ( I ,J ♦! I^va ,J*1 , KI-3*HT( I, J I >f\/ ( I , J ,K I -H7 ( I , J« 2 I * 


57 


CVIItJ«29Kn/(2*90Y) 

SB 


GO TO 2H 

59 

19 

COKUNUE 

50 


0lHUXr<4%HT<l4l ,J )9U II ♦ 1 , J, K 1 -3 '►HT ( 1 , J 1 C I , J > -NT ( I ♦ 2, Jl * 

61 


CU( I42f J, KU/( 2 .^ 0 X) 

5^ 


DIHVY=< 3*HT« I fJ )9VCI ,J»K )4H7 (l9J-2> *V(ItJ-2»K )-4*HTCl» J-U* 

63 


CVCI,J-l 9 Kn/( 2 »>^ 0 Y) 

5^ 


GO TO 24 

55 

16 

CONTINUE 

66 


01HUX = (3«HT( 1 »J {♦Ua t^,K MH T ( 1-2 9 JMU( 1-2 9 J i K 1-4 4rHT Cl -1 9 J 1 9 

67 


CUf I-l9J9KII/(29i^0X) 

6 S 


DlHVVriM •HTCI 9 J ♦ V CI 9 J ♦ 1 1 K ) -3 *HT f 1 9 J ) C 1 9 J fK ) -HTC I 9 2 » ♦ 

69 


CVCl 9 J« 29 K))/( 2 .« 0 Y) 

TO 


GO TO 24 

71 

17 

CONTINUE 

7Z 


D1HUX=C3*HTCI 9 J )»UC I , J ,K 1 ♦H T < 1 - 2 , J> ^U C I - 2 9 J »K I-4«HT « I - 1 , J) ♦ 

73 


CU<I-l 9 U 9 «n/( 2 .tDX) 

7M 


OlHVY-t 3*HT(I 9 J l^VCI 9 JfK I ♦H T C 1 9 J-2 1 *V ( I 9 J “2 9 K ) -4 <^HT C I 9 J*1 > * 

75 


CV(l 9 J-l 9 K)l/( 2 #« 0 YI 

76 


GO TO 24 

77 

24 

CONTINUE 

7B 

c • • 

••SIMPSON’S RULE IS USED FOR INTEGRATION 

79 


irCKtEQ^ 1«0R«K.E0.: ) GO TO 101 

80 


IF(K^EQ*2«0R»K«E0^4) CO TO 102 

81 


H 0 UH(l 9 J|r( CO lHUX«DlHVY)^nZi»C2^/3^) MHDUM (1 9 Jl 

82 


GO 70 103 

63 

101 

HOUHlIv J)=MO lHUX«01HVY)«0Z/3» )«H0UMCl9 Jl 

64 


60 TO 103 

35 

102 

HOUMCIf Jl=( (0 lHUX401HVy)«0Z*C4^/3.l MHOJHCI 9 Jl 

86 

103 

CONTINUE 

87 


60 CONTINUE 

88 


HTEC If Jl =HTO( I 9 Jl -HDUM (I ,U) ^2^0J 

B9 

60 

CONTINUE 

93 


RETURN 

91 


END 


90 


7.1.10 UVVEL 

This subroutine calculates the horizontal components of 
velocity, u and v, at t • At (-HCI.J.K) and G(I,J,K)), res- 
pectively) from u and v at t"0 (-U(I,J,K) and V(I,J,K))by 

using a forward differencing In time. 

Volume I details how the u and v momentu.Ti equations are solved. 
Note, this subroutine Is used only for the first time cycle. 

The general Inlet and outlet conditions are specified by reading 
In from data cards values of the parameters which set the location 
properly or the boundary. Modification, as mentioned earlier 
In the description of subroutine HEIGHT, may be easily Incorporated 
by the user. 

The spatial derivatives have been replaced by central diff- 
erencing in the interior of the domain and three point single 
sided differencing on the boundaries, as shown in Volume I. 

Again, MAR=6 and MAR=»8 boundary corners are tested as interior 
points . 
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C 

c********* *********** ***••«*»« 4>*»*»*«i«*********4t**»4i*#*******«*****ili»4M>**-^ 

C THIS SUBROUTINC CALCULATES THE HORIZONTAL V CLOCI TIES* U, V, AT EACH 

C X-V LOCATION AND DEPTH IN THE DOMAIN FOR THE FIRST TIME STEP USING 

C A rORUARO OIFFEPENCING SCHEME IN TIME 

C** *»**■»**»■•«»•*******•***««*«*»* ****«' ****»4>***4i»>t>*»««4««*<i**** 4=*** <■««' **<!•< 
SUBROUTINE UVVELflN *UN,KN ,U »V*H*6»0X*DV »OZ»W *TAUX*TAJY* 
COTtHTtHX*HYtETA,P*MAR*CItCC»CP»CH*CV»KH*KV»CR»RR (FFtHTD, R0» T* I 1 
CfI2* 13* 14*15* 16 *J1*J2*J3* J4 * J5 * J6 *H t 
REAL KH*KV 

DIMENSION UtIN* JN*KN)*V(IN* JN*KN)*H(IN*JN*KN)*G( IK'*JN*KN>* 
CHTIIN*JN>*HX(IN*JN) *HV (IN * JN ) *ET A ( IN *JN ) * 

CP(1N*JN*KN)*MAR (lN*JNI*HTD(IN*JN>*y (IN*JN*KNI 
C*R0( IN* JN*KNI *T (IN*JN*KN > 

KNX=KN-1 
DO 1C 1=1*IN 
DO 10 J=1*JN 

IF(MAR(I*J>.EO.OI GO TO 10 
IF(MAR(I*JI.E0.5> GO TO 10 
IF(HAR(I.J).E0*71 GO TO 10 
IFCMARd *JI .E0.9I GO TO 10 
IF(MAR( I * JI.EO. 10) GO TO 1C 
IF(MAR(I*J1.E0.3) GO TO 10 
IF(MAR(I*J).EC.4) GO TO 10 
DO e k:i*kni 

IF(MAR(I*JI.E0.6) GO TO 11 
lF(MAR(I*JI.EO.e> GO TO 11 
IF(MAR(I*J1.E0. Ill GO TO 11 
IF(MAR(1*J|.E0.1I GO TO 101 
IF(HAR(I*J).E0.2I GO TO 1C2 
11 CONTINUE 

ETAX = IETA(I*1,J )-ETA(I-l,J) )/ (2*0X1 
ETAY=IETA(I*J4l )-ETA(I*J-ll >/(2*0YI 
OlPXrCPC I-»l*J,K l-P(I-l,J ,K> )/ (2*0X1 

01HUUX = ( U(I*l*J ,K )*U'.I-»1 ,J*K)*HT (!♦ 1*J)-U (I-1*J*K)*U( I-l* J*KI* 
CHT(I-l*jn/(2*DX) 

D1HUVY=( U(I,J*1 *K )*V(I ,J*1*K)*HT(I* J«1)>U (I , J-1 , K )*V( I * J-1* K » * 
CHT(I*J-1)»/(2*0YI 
01UX = (U( I*1*J*K l‘'U(I-l*J*K) 1/ (2*0X1 
D2UX = (U( I«1*J«K t«U(I-I,J*K)<-2*U(I.J*Kn/(0X*0X} 

D1UY = (U( I*J«i*K I'UdiJ-ltK) )/(2*0V) 

02UY = (U( I*J*1*K )«U(I*J-l,K)-2*Uil*J*Kn/(0V*0Y) 

GO TO 100 
101 CONTINUE 

IF( J.EO. J1*AND. I.GE.Il.ANO. I.LE.I2) GO TO 10 
DHYr(3*HT(I*J »4HTII , J-2 > -‘••HT ( I * J-1 »)/(2*0Y » 

ETAXZ(ETA(I4 1,J )-FTA (I-l ,J» )/(2*0X » 

ETAY = (3*ETAI1 ,J )-»ETA(I , J-2) -4*ET Ad ,J-in/(2*0Y) 

D1PX = (P( I4 1,J,K )-P(l-l,J,KM/(2*0X) 

D1HUUXZ( U(14 1 ,J ,K l*U(I 4l , J,K »*HT (I* 1 * J I -U ( I -1 , J , K l*U( I-l * J* K I * 
CHT(I-l,Jn/(2*DX> 

DlhUVY=( 3*U(I ,J ,K )*V(I ,J ,K)*HT( I , J ) *U ( I , J -2 , K ) *V ( 1 , J-2, K I 
C*HT( I, J-2)-4*U( I, J-1 ,K )»V (1 , J-1 ,K J*HT (I, J-1 ) » /(2*DY) 

01bX = (U( I-*l, J ,K »-U(I-l ,J ,KI) / (2*0X1 
0 2UX = (U( I'1,J l«U(l -1 ,J ,ld -2*U ( 1 ,J ,IO ) / ( CX*OX ) 

D1UY=0.0 


93 


S7 


02UYS(UtX»J»K )«U(X,J-2|K t>2 *U(XfJ-l tKI)/(OV*OYl 

SO 


60 TO 100 

S9 

102 

CONTINUE 

SO 


XFU.E0.J2*AN0»I.Ce*X3.AN0.I.Le.I9| CO TO 10 

SI 


ETAX?(eT*a«l,JI>ET«(I-l tJ! )/ (2*0X1 

S2 


ETAY = (4*ETA(X tJ*n-3*ETA(X»J)-ETA(X |J« 2 n/( 2 * 0 V) 

S3 


0HYr(4*HT(X*J«ll>3*HT(XtJl*’HT IX(J«2 n/(2*£Y I 

s« 


OlFXSIPt X*1»J«K l•P(I-l,J,KI )/ (2*0X1 

ss 


01HUUXS(U(X«ltJtK l*U(X«lyJtKl*HT (X«lf JI'U(X-ltJ»Kl*U( X-l( J(KI* 

ss 


CHT(X-lt J)1/(2*0XI 

ST 


01HUVYS(9*U(X yJ«lyK}*V(Xf J«l,K}*HT(I,J«l)-3*U(Xt>itK)*V(Xt J,K1 

SO 


C*HT(XyJI -U(Xt 2tK)*V(l * J«2 yKI*HT (IyJ«2l l/( 2*0Y) 

S9 


01UX:(U( X«ltJfK )-U(I*l|JtK) 1/ (2*0X1 

73 


02UX=(U( X«ltJ*K 1«U(I*1|J tK) -2*U(ltJ,K}l/(DX«0Xl 

71 


OlUYsO.O 

72 


D2UY = (U( XyUyK l-»t'(I,J«2,K 1*2 *U ( I » J«1 yK 1 1 / ( 0Y*0V 1 

73 

100 

CONTINUE 

79 


RO(I yJyKl =1.0299 31-.00C320*T(XyJyKl-. 0000098* (T( I yJyXl **21 

75 


RRSROdy JyKl 

7S 


XF(K.E0.11 60 TO 70 

77 


01UWZ=(U(XyJyK* 11*U(I yJyK<*l 1«U(I yJyK>ll*M (I yJyK>lll/( 2*021 

78 


02UZ = (U( IyJyK«l l«U(IyJyK>ll*2*U(IyJyKll/(DZ«0Zl 

79 


60 TO 80 

00 

70 

DlUWZ=(9*U(Iy JyK«ll*W(I yJ,K4ll-3*U( XyJyKl*U (IyJyKl-U( Xy JyK«21* 

81 


CW(Iy JyK421)/(2.*0Zl 

02 


D2UZ = (2*U(IyJyK4l)4(TAUX/KV 1*HT (I yj l*2*DZ>2*U (I y JyKll /( 0Z*DZ1 

03 

80 

CONTINUE 

89 


UX=CI*(0 1HUUX40 lHUVY4HT(IyJl*01UbZl 

85 


UP=CP*HT(IyJl*(CTAX*GRl* (-1 . 1 

es 


UC=CC*HT(XyJl*Fr*V(IyJ,K 1 

87 


UH=CH*KH«(HX(Iy J1«01UX4£TAX*01UX«HT (I , J 1 *02UX 1 « CH*KHM HY ( I y Jl *0 lUY 

88 


C4ETAY*01UY4HT(I fJ)*02UYI 

89 


UV=CV*KV*02UZ/H TdyJl 

90 


H(IyJyKl= (('UI^l'P-UC^UH^UV l*OT«HT(lyJ)*U(IyJyKl }/HT3(Iy J1 

91 

8 

CONTINUE 

92 

10 

CONTINUE 

93 


1F(M.EQ.11 GO TO 7000 

99 


IFdNLET.EO.ll GO TO 7000 

95 


00 6000 K=lyKNl 

96 


SUH=0.0 

97 


00 6002 J=J3yJ9 

98 


SUH=SUH4H(I5*ly JyKl 

99 

6002 

CONTINUE 

130 


00 6003 J=J3yJ9 

131 


H(I64lyJyK)=SUM/(J4-J34ll 

132 

6003 

CONTINUE 

133 

6300 

CONTINUE 

139 

7000 

CONTINUE 

135 


00 970 I=lyIN 

13S 


DO 970 J=lyJN 

137 


00 960 f =1,KN1 

109 


IFd.EQ.IS.ANC.J.GE.JI.ANO. J.LE. J9. AND.M.GT.n CO TO 5002 

139 


IF(I.EQ.I6.AND.J.6E.J5.AND. J.LE. J6) CO TO 2CG5 

110 


60 TO 970 

111 

5C02 

H(I,J,KI=H(I4 1, JyKl 

112 


GO TO 970 

113 

2C05 

Hdy J,KI =Hd*l, JyKl 




u« 

US 

US 

117 

iia 

119 

iza 

121 

122 

123 

12<i 

125 

12S 

127 

12B 

129 

133 

131 

132 

133 
139 
135 
135 
137 
133 
139 
143 
141 
14? 

143 

144 

145 

146 

147 
14B 

149 

150 

151 

152 

153 

154 

155 

156 

157 

153 
159 

150 

151 

152 
163 

154 

155 
155 
157 

153 
159 
170 


960 CONTINUE 
970 CONTINUE 

00 30 I=ltIN 
00 30 J=1»JN 

ZF(MAR(I,J>.C0.0) 60 TO 30 
ZF(H*R(I,Ji.E0.5> 60 TO 30 
ZFIMARd (J> .E0.7I 60 TO 30 
ZF(NAR(I»J| .E0.9I 60 TO 30 
IFCHAR(l,Jt.E0.10l CO TO 30 
DO 7 K=1«KN1 

IF(HAR(I,J).E0.6> 60 TO 12 
ZFIHAR(I,J1.E0.(3) 60 TO 12 
XF(MARtI»JI.E0.11l 60 TO 12 
ZF(HAR(I»J).EO. 1) 60 TO 201 
IFIHAR(I*JI»E0.2> 60 TO 202 
IF(MARtI»J).E0.3l 60 TO 203 
XFIMAR(I,J).E0.4) 60 TO 204 
12 CONTINUE 

iTAX=tETAII«l,J )-ETA(I'l,J) )/(2*0X) 

ETAY=»ETA<ItJ«l >-ETAtI ,J-1) >/(2*0Y> 

01PY=<P< I,J*1,K »-P(I ,J-l ,KJ »/ (2*0Y» 

D1HUVX=(U(1«1*J »K )*V(I«1 t J*K)»HT (I« 1,J)-U (1*1 f JfK)*V( I*lf J*K)>k 
CHTtI-l,JJ)/«2*0X» 

01HVVY=( V(I»J«1 ,K )*V(I ,j4i,K|*HT (I, J4D-V (X ,J-1|K)4<V( I, J-1,K)* 
CHT(ItU-in/(2*0Y> 

OlVXrivt I4lt J»K l-VIl-lfJtKn/CZ't'DX) 

D2VX = (VI I«lfJtK )«V(I-l,J,K)>24>V(l,J,K)l/IOX*DXI 

01VY = (V(I|J« 1»K )-VII ,U-1 |K ) )/ (2*0YI 

02VY=(VI ItJ4l»K >«V(l,J-l,K)-2«V(I,J,Kn/(DV*0Y) 

60 TO 200 

201 CONTINUE 
60 TO 30 

202 CONTINUE 
60 TO 30 

203 CONTINUE 

IFil.EO. I5.AND. J.GE.J3.AND. J.LE. J4 ) 60 TO 30 
ETAX = (M*ETA(I41 ,J »-3»ETA (I, Jt-ET A(l42tUJ ) /( 2*0X I 
ETAY=IETA(I,J41 )-ETA(I,J-l J )/ (2*DY) 

0HX=<4«iHTIl4l ,J l-SfHTII , J )-HT { I ■» 2 i J M/(2*DX J 
OIPYSIPI I,J4 1,K 1-PCI ,J-l,K J 1/ t2*‘0Y) 

01HUVX=( 4«UCX«lf JfK 1*V CI41( J,K)«HT(l4l t Jl-3»UtI » JfKIAVC It J*X) 
C*HTIItJ) -UCl* 2t JtKI*VCI«2 fJ |K 1 «HT ( I 42 , J ) 1 / ( 2*0X 1 
01HVVY=C VCIt J4ltK )4<VCI |J4l,K I4HT (It J4l l-V (I t J-1|K)*V( I| J-ltK 1* 
CHT(ItJ-nW(2*0Y» 

D1VX=0.0 

02VXr(V(ItJtK l4V(l42tJtK 1-2 4V(l4ltJtKI)/( 0X00X1 
OIVYZCVC :tJ4l,K 1-VCltJ-lfKI 1/(200 Y1 
D2VY = (V( I,J4ltK |4V(ItJ-ltKl-24V ( ItJtKII/ (DYOOYI 
GO TO 200 

204 CONTINUE 

IFCI.EQ.I6.AND. J.CE.J5.AN0.J.LE.J51 GO TO 30 
ETAX = ( 30ETACI ,J 14CTACI-2 , J1 "40ET A( 1 -1 fJl 1 / C 2O0X 1 
ETAYr«ETAII,J4i )-ETA(I tJ-11 5/ ( 2 oOY) 

0HX=(3oHT(I ,J 14HT «I -2,J1 -40MT (I-l ,J 1 1/C20DX 1 
01PY = (P( l,J4i,K )-PCI,,>-l,l<) )/(20DYl 

DIHUVX=C 3 4U(I,J,K )«V(I ,J ,K } oHT C I t J > «U ( I -2 tJ » K )0 V ( I J, K 1 
COHTC 1-2, J1-40UC J^Vd -1 ,J,K lOHTCI-l ,J 1 l/JZOOX) 
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in 

172 

173 
17« 
175 
17S 
177 
179 
179 
183 
191 
182 
19 3 
10« 
195 
IBS 
187 
19 9 
199 
193 

191 

192 

193 
19<l 
195 
19S 
197 
199 
199 
233 

231 

232 

233 
239 
235 
23S 
237 
203 
239 
210 
211 
212 

213 

214 

215 
21 S 
217 
213 

219 

220 
221 
222 
223 


CHT(X»J-m/(2*OV) 

oivkso.o 

D2VX:tV( I,J*K l«V(X-2»J*K 1-2 9V ( I*1 « J (K I W I OX *0X ) 

DIVVSCVI I,J«lfK )«VtI»J-t»KI l/(2*0V) 

02VYS(V( I,J«1»K )«V(XtJ<'l(K)>2«V(ItJrKI)/fDY*0YI 
2C0 CONTINUE 

RO(ItJ»K ):U029431-.COC020«T(ItJ,K)-.00000489mif J,K79*2I 
RRSROUt JfKI 
inx.eo*!) GO TO 9C 

01VUZs(VCXtJ(K« l)<*wa»J«K«l l-V(X|J»K«l l*W (X,J»K'1I)/(2*0ZI 
02V2stV( X,J,K«1 }«VII ,J»K«1)-2*V(I ,J ,Kn/(0Z*0Z) 

GO TO 95 

90 D1VWZ:|4 WV(I,J»K«1)*WIZ»J,K«1)-3«V(I,J,K (XtJf Kl'Vf I« JtK«2>* 
CWIItJtK«2))/( 2.O0ZI 

02V7=(2*V(I»JtK-«n4(TAUY/KV I9HT II tJ )*2 *0Z -2* V (I t J tK II /1 02*021 

95 CONTINUE 

6 0 VI=CI*ID 1HUVX4D 1HVVY«HTI I »J )*01Vy7) 

VP=CP*HTII,J)*IETAY*GRl*l-l. I 
VC=CC*HT(I,JM<FP*U(I »J,K ) 

VHrCH*KH*(HX( I, J)*01VX«ETAX«>01VX«HT ( 1 » J ) *02 VX ) *C H*KH* I HY ( 1 1 J ) *0 1 VY 
C*ETAY*01VY«HTII ,J>*D2VYI 
VV=CV*KV*D2VZ/HT(I»J) 

Gilt JtKI = 1 1 -VI ♦VP4VC4VH4VV 1*DT4HTI 1,J)*W (I ,J»K1 >/HT3l It Jl 

7 CONTINUE 

30 CONTINUE 

IF(N.EO.l) GO TO 700 
IFIXNLET .GT.l ) GO TO 700 
DO 600 KriyKNl 
SUM=0.0 

DO 602 I=IlfI2 
SUHrsUM«6IItJ 1-ltK) 

602 CONTINUE 

00 603 I=IltI2 

Gilt Jl~ltK)=SUM/II2-Il«n 

603 CONTINUE 

600 CONTINUE 

700 CONTINUE 

00 97 in, IN 
00 97 J=1,JN 
00 96 K=1,KN1 

IFIJ.EO. J2.AN0.I.GE.I3.AN0.I.LE.14) GO TO 205 
IF<J*£Q.J1.AND.I«6E.11.ANU.I.L£.I2. AND.H.GT.l) GO TO 502 
GO TO 97 
2C5 CONTINUE 

Gilt JtK)rca,J4 1,K> 

GO TO 97 
502 CONTINUE 

ClI, JtK)r6(I,J-l,K) 

96 CONTINUE 

97 CONTINUE 
RETURN 
END 



’ t- ■ ■ ■ ' 
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7.1.11 UWELN 

This subroutine calculates the horizontal components of 
velocity u and v at time level n+l (>D(I,J,K) and E(I,J,K) 
respectively) from u and v at time level n (>>11(1, J.K) and 
Gd.J.K) and u and v at time level n-1 (?U(I,J,K) and V(I,J.K)) 
by using central differencing In time. The numerical scheme used 
for solving these equations Is given In Volume Again, the 

general Inlet and outlet conditions are specified and may be 
modified by the user. Note, DuFort-Frankel differencing is used 
for the vertical momentum diffusion term. The spatial derivatives 
have been differenced as shown in Volume I, 
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1 

2 

3 

« 

5 

6 
7 
• 
9 

10 

11 

1? 

13 

It 

15 

15 

17 

18 
19 
23 
21 
22 

23 

24 

25 
25 
27 
23 
29 
33 

31 

32 

33 

34 

35 
35 

37 

38 

39 

40 

41 

42 

43 

44 

45 
45 
47 
4 3 
4 9 
5D 

51 

52 

53 

54 

55 
55 


C 

C***««r**»>l***«*M*««4M»«*****»9******»««*»**9*****99«**********9********9 

C THIS SUOROUTIK'Z CtLCULATCS THE H0t7I70NT«L VELOCITIES* UfV* AT EACH 

C X-V LOCATION AND DEPTH IN THE OOHAIN FOR THE SECOND IINE STEP ANO 

C THEREAFTER USING A CENTRAL DIFFERENCING SCHEHE IN TINE 

C*****9**«4*»* ******* *«M***4>»*»**»***<»*«**9****»***9* **«»«*»«*«* •* 

SUBROUTINE UV VEIN <1N »JN ,KN * U,V *H *C * D»E * DX »0Y , 02 » U * T AJX* T AUY * OT * 
CHT,HTO,HTE*HX»HY*fTAtPfHAR,KM«KVtGR*RR»FF tCPtCCtCI *Cr4 * C V , RO » T *1 1 
C»I2* 13*I4»IS,I6tJltU2tU3(J4 , J5 * J6 *H ) 

REAL KH y K V 

OIPENSION UdNyJNyKNlyVIINyJNyKNlyHtlNyJNyKNIyGt INyJNyKNI y 
COdNyJNyKNIyE (IMyJNyKNIyHTI INyJNlyHTOdNyJNlyHTEdNyJNI yHX(INyJN) y 
CHYdNy JN I yLTA (I r;y JN) yP(INyjr4yKN» yHARdNy JN) yWdNyJNyKN) 
CyROdNyJNyKN) yT (INyJNyKNI 
KN1=KN>1 
00 lU iriyZN 
00 10 JZlyJN 

IFiMARdyJI .EO.n GO TO 10 
IF(MARdyJ) .EO.'j) SO TO 10 
IFtMARd yJI .EO. 7) CO TO 10 
ZFIHARdyJ) .EC.9I 50 TO 10 
IFtMARd yJI .EC. 10) CO TO 1C 
IFCNARdyJI.CO. 3) GO TO 10 
IFIMARdyJI.EC.t) (iO TO 10 
DO 8 K:1 yKNl 

IFtHARd yj) .EO. 5) GO TO II 
IF(MARd yJ) .EO.C) GO TO 11 
IF(NAR( I y J) .EC. 11 1 CO TO 11 
IF(MARd y J) .EO. n GO TO 101 
IF(f*AR( I *J) .EC.2) GO TO 102 
11 CONTINUE 

ETAX = CETA(I*1,J l-ETA (1-lyJ) )/ (2*0X1 
DHX=CHT( I*lyJ )-HT(I-lyjn/ (2*0X1 
OHY=(HT( IyJ*l »-HT(I yJ-1) )/(?*DY) 

D1PX = (P( I*1,J,K )-Pd-l,J,K) )/(2*DX) 

01HUUX=(H(I4 1 yJ.K )*H(I«1 yJyK )*HTO(I ♦lyJ)>Hd«lyJyK)« 
CHd-lyJyK)*HTrd-lyJ 11/ (2*0X1 
01HUVY=(H(I,J«lyK >*G (I yJ«lyK)*HTO(I yJ^U-HdyJ-1 yKI* 

CGdy J-ly K)*HTDdyJ-l)l/(2*0Y) 

01UX = (U( I«lyJyK )-Ud-lyJyK) )/(2*DX) 

02UX = (U( I«lyJyK )«Ud-lyJyK)-2*UdyJyK)l/(0X*0X) 

DIUY=(U( »-Ud ,J-1,K) )/(2*0Y) 

« 02UVr(U( IyJ«l yK )«UdyJ-lyK)-2*UdyJyK))/(0Y*CY) 

GO TO 100 
ICl CONTINUE 

IF(J.EO. J1.AN0.I.PE.I1.AN0.I.LE.I2) GO TO 1C 
ETAXZ(£TA(I*1,J l-ETA (I-1,J) )/(2*UX) 

DHXr (HT( Id, J )-HT(l-l,J» ) / (2*0X1 

OHY=( 3*HT(I ,J l + NTd ,J-2) -‘(■t-HT (1 n/l2*DY I 

DlPX = (P(I*l,J,K)-P(I-l,J,Kn/(I’«L)X» 

D 1HUUX=( H(l4 1 ,J ,K )*H (T ♦ 1 , J,K )*HTU (I ♦l,J»-Hd-l,J,K)* 
CH(I-l,J,KI*HTO(I-l,jn/(l'*OXI 
D IHU VY= ( I«H( I ,J ,K l»6 (I t J ,K I ■>HTD( 1 ,J ) «Hd , J-2 yKI*G( I , J-2, K I 
C*HTP (I,.J-2) -<l (1 ,J-1 ,K I *G( T ,.J-1 ,K ) ^HT D ( 1,J-1))/(2»UVI 

0 1UXZ(U( IV 1, J ,K )-U( I-l , J,K : )/ (2*CX I 
D2UX=(U( ivlyd ,K )«U( I-l ,J ,K I -2*U ( I yj ,KI )/(DX*DX) 


> 4 »* 


101 


57 


OIUYSO.O 

59 


02UY:IU( I»J*K)«U(I»J-2(K l-2*U(Z*J-l«KII/IOY*OY> 

59 


60 TO 100 

50 

102 

CONTINUE 

bl 


ir(J.EO.J2*ANO.X*6E.X3.ANO.ZtLE.Xb> CO TO 10 

b2 


ETAX:(ETAfI«l»J»-CTA(X*l|J) )/( 2*0X1 

bS 


ONXSIHTI I«ltJ)-HTfl>ltJ> >/< 2*0X1 

Sb 


OHYr|b*HT(l,J*l)-3*HT(Z»Ui-HT(I»J«2 n/(2*0Y » 

bS 


D1PX:(P< I«l(UtK l•P(X•l(J»K) )/(2*UX) 

bb 


01HUUX:(N(X«1,J(K l»H (X«l t J»K}*HTD(l«l»JI>Hf X-ltU»K>* 

b7 


CH(1-1»J»KI*HT0(I-1,J)I/(2*0X) 

69 


OlHUVYSf b*H(I »J«1,K )*G(I yJ« lyK>*HTD(X*J«n>3*H(X li J(K ) 

69 


C*HTD(IyJI*HIZ yJ«2yK i*G(I yJ«2yK >*HT0llyJ*2))/(2*0YI 

73 


01UX:(U( I«lyJyK )-U( I •! yJyK ))/ 12*0X1 

71 


02UX = (U( I«lyJyK )«U(I-lyJyK) >2*011 yJyK)l/IOX*OXX 

72 


01UY=0.0 

73 


02UYSIU< ZyUyK l*UfIyJ«2yK )-2*U(XyJ«l yKI l/(OY*OY) 

7b 

100 

CONTINUE 

75 


RO(IyJyKlsl.O29b31-.00002O*T(ZyJyK>-.OO0OOb8*(T( XyJ»X)**2) 

76 


RR=RO(Iy JyKI 

77 


XFtK.EQ.ll GO TO 70 

73 


01UU2=tHIIyJyt<« l)*b flyJyK*! l-HlIyJyK'l 1 *U (X yU yK • 1 1 1 / f 2*0 2 1 

79 


UZ=H(ZyJyK«ll«H (IyJyK>ll*Uf lyJyKl 

80 


60 TO 80 

81 

70 

01UU2 = fb»HayJyK«l)*U(IyJyK«l )-3*H(IyJyKl*lMIyJyKI>Hf ly JyK*21* 

82 


CU( lyJyK* 2)l/( 2. *02) 

83 


UZ:2*H(IyJ»K« 1) «(rAUX/KV)*H70(XyJ)*2*0Z>U(I yUyKI 

8b 

6C 

CONTINUE 

85 


UI=CI*IO IHUUX40 lHUVY«HT0IIy J)*01UNZ 1 

8b 


UP=CP*HTOayJ)*(ETAX*GRl*{-l. 1 

97 


UC=CC*HrD(IyJ)*FF *6(IyJyK 1 

88 


UH=CH*KH*OHX*D 1UX*H T( I , J >*02UX ) ♦CH*KH'MDHY*01UY «MT(Z yJl *02UY ) 

89 


t](IyJyKl = (( (>UI ♦UP-UC*UH 1*2 *0T*» 1CV*KV *UZ 1/ 1 0Z*0Z*HT3 1 1 y J 1) 1*2* 

90 


C0r«U(IyJyKI*HT( I,JH/HTE( ly J) 1/ (l.« (Ctf*KV/(CZ*DZvHTO( ly J) 11*2*07 

91 


C/HTE(IyJll 

92 

B 

CONTINUE 

93 

10 

CONTINUE 

9b 


IFIN.EO.ll GO TO 700C 

95 


IFdNLET.EO.l 1 CO TO 7000 

95 


00 6000 K=lyKNl 

97 


SUH=U.O 

93 


00 6C02 J=J3yJ4 

99 


SUH=SUM«0(I5« ly JyKl 

130 

6002 

CONTINUE 

101 


00 6C03 J=J3yJb 

102 


0(I5*lyJyKirsUM/(J4-J3*l 1 

103 

6003 

CONTINUE 

10b 

6000 

CONTINUE 

135 

7C00 

CONTINUE 

106 


DO 970 I=lyIN 

137 


00 970 

139 


00 960 K=lyKM 

139 


IF(I.EG.I5.AND.J.GE.J3.Af.n.J.LE.J4.Ar,D.*f.GT.ll GO TO 5300 

113 


IF( I.EO. le.AND. J.CE .J5 .AK'D. J.LE. J6 1 60 TO 2005 

111 


GO TO 97C 

112 

5 002 

0( I, JyKl rO<I« XyJ.K) 

113 


GO TO 970 
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ii« 

ns 

ns 

nr 

ns 

119 

129 

121 

122 

12S 

129 

125 

126 
127 
12S 
129 
139 

131 

132 

133 
139 
135 
13b 
137 
1S9 
139 
199 

m 

192 

193 
199 

195 

196 

197 
195 
199 
159 

151 

152 

153 
159 

155 

156 

157 
153 
159 

150 

151 

152 

153 
159 
155 
155 
157 
153 
159 
170 


$ 


2C05 0(Z«JtKt:OII-lt«ltM 
96C CONTINUE 
970 CONTINUE 

00 30 isltXN 
00 SO J=ltUN 

ir(NAR(Zt*ll*EC.C> GO TO 30 
ZFfNAR(Zt«ll*EC.5> 60 TO 30 
ZFtMAR(Z»J)*E0.7> 60 TO 30 
ZFlMARtZf J).CC«9) GO TO 30 
ZnHARtZfJl.EO.lOl 60 TO 30 
00 7 KSltKNl 

Zr(HAR(Z»Jl.C0.6l 60 TO 12 
ZFCMARIZtJI.EO.ei 60 TO 12 
IF(HAR4It«n.E0«lll 60 TO 12 
IF(MAR(I,J).eO. n 60 TO 2C1 
IF(HAR1I,J).CC.2I 30 TO 202 
ZFIMARII ,J)*e0.3) GO TO 203 
ZFCMARIZ tJ>.E0.9) GO TO 209 
12 CONTINUE 

etay:(etaii»j 4 1 1-ETni ,j-m/c2*oyi 
OHXSiHTt I«1»UI-HTII-1»J> l/(290X> 

OHVsiHTt XiJ«l t^HTII |U«1) )/(290V) 

OlPVrCPI >-P(l ,J-1,K> )/(2*0VI 

OlHUVXsf H(Z« 1 tJ tK >*G (Z« 1 tJ(K )*HTOn«l»JI>HI l-l|J»K)*6f I-I* JfK I* 
CHT0(Z-1» JI)/(?*DX> 

01NVVYstC(l(J«l*K 1*6(1 ,J«ltK)*HTD(l tU«ll*6( ZtJ’l »K)* 

CGUf J-lf K)*HTD( 1»J-1 ) ) / (2*0YI 
01VX=(V( Z«l»JtK I V(I-1»J,KI)/(2*0XI 
D2VX=(V( I«lyJ ,K )«V(I-1 ,J ,K) -2*V ( I , J «KI } / ( OX *CX ) 

01VV=(VI 1,J*1,K »-V(I,J-l ,K»J/(r*OV» 

02VY=(V( I«J«1»K l*VII .J«l tK)-2*V (I*J»KM/(OY*OY) 

GO TO 200 
2C1 CONTINUE 
60 TO 30 

202 CONTINUE 
GO TO 30 

203 CONTINUE 

IFII.EO. lE.ANH. j. 6C.J3*AND. J.LE. J91 GO TO 30 
ETAY=(ETA(I,J*l)-ETA(I ,J-1) 1/ (2*DV) 

0HX=(9*HT(I*l,J)-3*Hni, J)-HTII*2»JII/{2*0X I 
DHY=(HTI I,J*ll-HT(I,J-l» »/(2*DY) 

01PY=(P( I,J*1,K l-P<I ,J-1 ,l<n / (2*DYI 

OlHUVXSt 9*H(1 «1 vJ«K )*6 (7*1 , JfK >»HTD (l«lf JI-3*H(I »J»K>*3( If J«K ) 
C*HT0(Z,J)>H(I«2 fJfK )*G(I«2,JfK )*HTO (1*2 fJli / (2*0X1 
01MVVY=( G(I*J«1 fK )*3(1 fJ*lfK l*HTD(I fJ«lt«6(IfJ-l (K)* 
CGII,J-l,K»*HTD(I,J-in/«;:*DY) 

01VX=0.Q 

0 2VX = 4W( IfJiK »♦ V(I*2,J,K »-2*V<I*l,J,Kn/(DX*DX» 

D1VY = (V( IfJ* 1,K l-V(I f J-1 ,K 1 )/ (2*0Y1 

02WY = IV( If J* 1,K >*V«1 ,J-1 fKJ-:*V«lfJfKM/4 0Y*DY» 

CO TO 200 
209 CONTINUE 

ird.EQ. I6.ANC.J.GC.J:.AN0. J.LE. J6) GO TO 30 
ETAY=(ETAIIf J*1 )-ETA (I ,J-1) »/ (2*DYI 
DHX=l3*HTa t J )*HT(I -2f Jl -9*HT (1-1 ,J »>/<2*DX 1 
0»IY=IHT( IfJ*n-MT(IfJ-lU/(r*DY» 
DlPYr(p|IfJ*l,KI-P(IfJ-lfK))/(2*DYJ 


INTCN-T^OKWiW. 


105 


111 

172 

175 

17* 

175 

176 

177 
171 
179 
119 
111 
112 
113 
119 
IIS 
111 
117 
HI 
119 
193 

191 

192 

193 

194 

195 

196 

197 
19S 
199 
230 
201 
292 
233 
294 

235 

236 

237 
236 
299 
213 
211 
212 

213 

214 

215 

216 

217 

218 
219 
22 0 
221 
222 

223 

224 

225 


OlNUVXri 3 *HI 1 tJ tN H 6 IItJ*K 19 HTD(l»JI«H(I' 2 f J»K»*GII* 2 »«>»KI 
C*HTO(X- 2 tJ>- 4 «N II -l,JrN )»CI I-l ( J , K I *HT Of I - 1 1 / 1 2 *OX I 
01 HVVVSl 6 ll(J«ltK>*CII tJ«ltK MHTOIl tJ«l 1 -GI I (J -1 
CCIlt«f<-ltKHHT 0 IXtJ*in/l 240 YI 
OlVXSO.O 

02 VX:lVI ItJtN )«VII- 2 »J(KI* 2 *Vli-liJ»KII/ICX*OX) 

OIVVSIVI X,J« 1 »K l•V(l ,J-l,Kn/l 2 * 0 Yt 

D 2 VYSIVI X»J« 1 »K l•VtI,J-l»K)- 2 *V(I,J»Ktl/IOY*OY} 

2 CO CONTXNUE 

ROIl«JtNHl*O 29 ii 31 -*OOO 02 O*riXtJt*<l-.C 0 OCO 4 e«ITIt»U»K)** 2 l 
RRSROIXi J«KI 
XPlK.CQ.ll 60 TO 90 

01 VWZ:i 3 II»JtK« l)*wn,JtK«i I-GII »J|K- 1 I*W IX (J»K- 1 M/I 2 * 02 ) 

VZsCIXf J»K«1)«6 IXt^IfK-ll-VI IiJfK I 
60 TO 95 

90 01VUZ=l4*CIXtJt*<«H*<^(If JtK«l t-3»GI Xt«l*KI*U (X»J»NI>GI Xt J|K«2|* 

CWtI,J,K« zn/l 2.90Z) 

VZ: 2 *GIItJtK« IMITAUY/KV l*HT 0 IZ|J>» 2902 >VII»JtK) 

95 CONTINUE 

VX=CX*I 0 1 HUVX «0 1 HVVY«HTDII, J)* 01 VWZ ) 

VPSCP*HTC{ZtJ 1 * IETAV* 6 RM'I> 1 *) 

VC=CC 9 HTDII»JI*rP«HII,J,K I 

VH=CH*KH*I0HX*0 1 VX«HT(I t JH02VX ) «CH»KH*IDHV *01VY «HTIXf J) 90 2VY ) 
CIZfitKlslll>VZ«VP«VC«VH)*2 90T«l ICV fKVAVZ 1/ I0Z«0Z*HT3II» jm*2*OT 
C«ViXtJtK )*HTIXt Jn/HTEII iJ) 1/ 1 1 . « <C V9KV / 1 OZ *0Z«<H TO I Xt JM i <*290 T 
C/HTE IXtJn 
7 CONTINUE 

30 CONTXNUE 

XFlH.EQ.l) GO TO 700 
XFIINLCT.GT*! ) CO TO 700 
00 600 K=ltKNl 
SUHso.O 

00 602 I:I 1 ,I 2 
KUM:SUH«ElltUl*ltK) 

602 continue 

00 603 I=I 1»12 
ElI,Jl-l,Kl=SUM/II 2 -Ilil) 

603 CONTXNUE 

600 CONTINUE 

700 CONTINUE 

DO 97 1 = 1 , IN 
DO 97 J= 1 ,JN 
00 96 KSltKNl 

IFI J.EC. J2,AN0. I.Ge*I3,AN0* I#LE*I4 ) Go To 20 S 
IFIJ.EQ. Jl.AND.I.GE.Il.AND.I.LE.IZ.AND.M.GT.n GO TO 502 
GO TO 97 
205 CONTINUE 

ClI,J,K|=C(I,J«l,Kl 
60 TO 97 
502 CONTINUE 

EII,J,K|=E|I,J-1,K) 

96 continue 

97 CONTINUE 
RETURN 
END 
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7.1.12 WVEL 

This subroucine calculaces the equivalent vertical vel- 
ocity 0(1. J.K) In the oi, 3. o ("X.y.a) coordinate system at each 
X. y location and depth a In the domain 0(I,J,K) at t ■ At 
(■W(I,J,K)) Is calculated from u, v, and H at t ■ At (•H(I,J,K)) 
G(I,J,K), HTD(I.J.K)) as shown in Volume I. Thereafter, 0(1, J,K) 
at time level n+l Is calculated from u,v and H at time level 
n+1 ("D(I,J,K)j E(I,J,K), HTE(I,J)) Simpson's Rule Is used for 
performing the Integration. 
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1 

2 

C 

C*0* 



s 

c 

THIS SUDROUTINC CALCULATES THE EOUIVALENT VERTICAL VELOCITY 

IN TMC 

% 

c 

COORDINATE SVSTEh AT EACH X«Y LOCATION AND DEPTH IN THE OOHAIN 

s 

c 

AT EACH TIKE STEP 


6 

7 

t 9 V 

SUBROUTINE NVEL (1K,JN,KN ,U,V,5,HT tOX ,DV *07 , M AP ,K , 1 1 , 1 2» J3t J9 ) 

s 


OIPXNSION U(IN,JN«nN |,VI1N»UN»KN ).HT|IN,JNI (ri (IN , JN «K Nl « HAR ( I N, JN) 

f 


KNirHN-1 


19 


00 lo U‘ 1»1N 


It 


DO 10 JsltJN 


17 


DUH79. 


15 


DO 9 K-1»KS 


1« 


IF(HAR(I ,J» tEO.ri 00 TO lO 


15 


IRHARdf Jl.CO* 11) 60 TO 11 


U 


iriMARII »J) fEC.6) Go TO 11 


17 


1F(PAR(X t j) .ec.()) 60 TO 11 


15 


IF(HAR(l,J»,EC.l,AND.I.GC.Il.AN0.I.LC.12tAN0.H*GT.l> 60 TO 

12 

19 


inHAR(Z»J).£0.3.AND.J.Ce*J3»AN0.J*LE.J9.AND.M,CT.il 60 TO 

12 

79 


lF(r(A<’(Ifi>*tT.ll) CO TO 10 


71 

11 

OlHUXr(HT(I«i,J l«U(I«:,J»KI-HT(I«l, J>*U(I-1 ,J|KI l/(2.«<0XI 


77 


01HVY=«HT(I ,J^1 |*V<1 ,J4l ,K ,J-n*V(I ) / ( 2.9DY) 


73 


GO TO 24 


79 

12 

CONTINUE 


75 


CALL UVEt.l(I(J,K,IN»JN,KN»UtV*HT,OX »0V »HAR« 0 IHUX ,01HVV) 


75 

29 

Continue 


77 


IFIK.EO. l.On.K.eO.5) GO TO 27 


75 


1F|K.E0.2*OR.K. ec.it) GO TO ?6 


79 


0uH=DUH«0z*(2./3. |*|0IIIUX«01HVY)/HT II,J) 


39 


GO TO 9 


31 

27 

DUK=DUN* (02/3 .) <M0IHUX4D1HVY)/HT (It j) 


32 


60 TO 9 


33 

28 

DUH=0UM«0Z*(9 ./ 7.)»(DlHUX«0lHVY)/HT (I,J) 


39 

9 

continue 


35 


WUO=Ot 


35 


DO e K=2fKN 


37 


IF(HAR(It Jt.eO.Cl GO TO 13 


39 


lF(HAR(I»J).CO.m 60 TO 111 


39 


IF(MAR( It J) .eo. 6) 60 7C 111 


99 


IF(MAR(ItU).£C.S) 60 TC m 


91 


lF(HAP(ItJ)*£0.1.AN0.I.C£.lI.AND.I.LE.I?.AND.M.GT.l) GO TO 

112 

92 


1F(KaR( I t J) «£C. 7. AN0.J.6E. J3.AN0. J.LE.J9 . ANO.H.G T.l > 60 TO 

112 

93 


IF(KAR(ItU).LT. U) GO TO 10 


99 

in 

DIhUX=(hT(14 1 ,J )•U(I4ltJtK)-HT(I-ltJ)•U(I-l tJ,K) l/(2.*0X» 


95 


01HUX1=( HT(1* 1, J>-»U(14 1, J,K-1»-HT (1 -1, J)*U( I - 1 , J ,K -l> ) / ( 2 .*0X ) 

95 


DlHvY=(HT(Itd*l )<!'V(I,J4l,K)-HT(I (I,J-l,K) »/(2,*DYI 


97 


D 1IIVY1=( HTr I, J4 1)*V (I t J*1 tK-1 )-HT (I tU-1 >*V( 1 1 J" 1 i H - 1 ) » / ( 2 . *0 Y I 

93 


GO TO 200 


99 

112 

continue: 


53 


CALL WVEL2(I,J,KtIN,JN,KN,U,V,HTtUX ,LY , K A R t 3 IHUX tClHVYt 


51 


CQ1HUX1,D IHVYl } 


52 

ZOO 

CONTINUE 


53 


ir(H.E:u.z) 00 to ici 


59 


IF(K.CC.9» 60 TO 101 


55 


lr(K,Et.3» Go TO lo? 


55 

102 

wuo=iiUij*oz*( < r, /.?.»A(r jhuxi*dihvyi 1 ♦d./s .)-^ (oiHux^oinwY) »/ 
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57 


CHTC1»J» 

SB 


GO TO 300 

59 

ini 

WUD = WU0*02*(I l,/3,|*(DlHUXl«DlHVVn ♦J2,/3.J*J0lHUX40lHVVII/ 

S3 


CHT(ltJ) 

51 


GO TO 300 

6? 

3C0 

CONTINUE 

ol 


WCI»J»K) =-UU0 40UH«>(K>1 )*OZ 

64 

8 

continue 

65 

10 

CONTINUE 

66 


RETURN 

67 


END 








preceding page dunk not filmed 


110 


7.1.13 WEU 

This subroutine calculates the differential terms, in diff- 
erenced form, in the definite integral for the equivalent 
vertical velocity, at each time step, from u, v, H at t**At, 
and at the time level n+1, thereafter. 
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IM 

15 

IS 

IT 

IS 

19 
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C****»»********9«*«****«*«*»****»**************«******************9*******v 

C This subroutine calculates the differential (OIFFERENCES} terms in the 
c definite integral for the equivalent vertical velocity at 

C EACH TIME STEF 

C9*9******99«>A«***«*«<****«**********«<»»***»****»*****4>«*«*«***»*****»**«***< 

SUBROUTINE WVEL1(I,J,K ,IN,JN ,KN,U,V ,HT,DX »0Y ,MAR ,01HUX,0 IHVY) 

dimension UtINf JN»KN)tV(IN( JN,KN>»HT i1N,JN) .MAR(IM,JVI 

IF(MAR(ItJ).E0.3) GO TO 12 

IFfHARd ,Ji .EO. SI GO TO I9 

IFIMAR' i,J) .E0.2I GO TO 13 

IF(MAR:i»JI.E0.1l Go To ?0 

1FCMARII,J|.E0.4I go to 14 

IF(HAR(I,J).E0t7) GO TO 1$ 

IFlMARd tJ).E0*9) GO TO 16 
IFIHARdtJI.EO.lOl Go TO 17 

12 CONTINUE 

01HUX:(4 *HT(I«1 »J I*U(I«1 tJfK )-3<liHTdtU)*U (I tU»K> oHTII«2» Jl« 
CUII«2,JrKn/(2.«0X) 

DlHVY=(HT(lf J«ll>«>VdtJ«l»K)-HT( I»J-l)*Vd ,J*lfKI l/d.oOYl 
60 TO 24 

14 CONTINUE 

OlHUX:t3«HT(I yJ )9U(I » J ,K } «HT (1-2 *UH>U (1-2 f K )>4*HT(1 *1 1 J1 * 
CU(I-l»J,Kll/( 2.«0X> 

01HVY=(HT(ItJ«l H>V( I,j«1,k)>HT( I»J-ll*vd tU'ltKn/(2.9DY) 

60 TO 24 

13 CONTINUE 

IF(J.E0.1.AND.I .GE.31.AN0.ItLE.33l GO TO 31 
OlHUX:(HT(I«ltJ )*U(I«liJ|K)-HT (I>1,J|*U(I-1 ,J,KI 1 /(2,«>0XI 
01HVY=(4*HT(I ,J^1 )*>V(I ,j4l,K)-3*HT( I,J1*V (I,J,K) -HT(It J*2I* 
CV(I,J«2,K)}/(2,*0YI 
60 TO 24 

31 CONTINUE 

01HUX=(HT(I«1 ,J l«U(l4l ,J tK) -HT (1-1, J)*U(I-1 , J ,K I 1 / ( 2. * jX I 
01HVY=(4*HT(I »J ♦! ItVCl ,J*1,X )-3*HT(I,J)*V (I ,J,K1 -HT(I,J» 2l* 
CV(I,U«2,Kn/(2.*OY) 

GO TO 24 
20 CONTINUE 

IF(J.C0.11.ANC.I.GE.7.AN0.I.lE.le,l 60 TO 32 
01HUX=(HT(I*1,J )*U(I-*1,J,K) -HT (I-l, J)*U(I-1 ,J,K1 )/(2.«'0XI 
01HVY=(3*HT(1 ,J)*V(I,J,K )4HTII,J-2)«t-V(I,J-2,KI-4*HT(I,J-l>* 
CVd,J-l,K))/( 2.40Y) 

GO TO 24 

32 CONTINUE 

01HUXs(HT(l4l ,J )>KU(I«1 ,J ,K)-HT(I-1, J)*U(I-1 ,J,KI l/(2.4>0X) 

01HVY=(3*HT(I *J >»V(1 ,J,K )4HTiI,J- 2) *VIItJ-2,KI-4*HT(I,J-ll* 

CV(I,J-1, 101/(2. *OYi 
60 TO 24 

15 CONTINUE 

DlHUXr(4 ♦HTII ♦! ,J )*U (I 4 1 , J, K I -3>CHT ( I ,J» *U (I ,J,K) -HT(I4 2, Jl* 
CU(l«2,J,K))/( 2.*0X) 

01HVYr(4*HT(I , J ♦ 1 1* V ( I , J« 1 , K I -3*HT I I , J1 * v ( I , J ,K » -HT( I , J* 2 1 ♦ 
tVd,j42,Kn/( 2,*DY» 

60 TO 24 
J9 CONTINUE 

DlHUXr(M»HT(I+l ,J )»ll (I •* 1 , J, K ) -3*HT ( 1 , J » *U (I,J,Kl-HT(l42,J)* 
CU(I42,J,K)»/(?.*0X> 


113 


57 


01HVY:(3*HT(I ,J )*V(I tJ tX ) «H7 1 X , J*2 1 «V U tJ'2 f Kl •4*HTf I » J-1 ) * 

SB 


CV( Xi J*l« KU/{ 2. *OYl 

59 


GO TO 24 

S3 

16 

CONTINUE 

51 


01HUX:|:»HTIX ,J l*U(I,JtK I«HT(X-2,J1 *U(I-2 ».* tK 1*4 *MT (I -If J) «> 

s; 


CUIl-lfJ*K)t/(3.»OX) 

63 


01HVY = (4*HT(I fJ«n»V(If J«1,K)-3»HTC 1 Jl^VII f JfK) -HTIX f J«2)« 

64 


CVtXf J«2fK)l/(2.*OY) 

65 


GO TO 24 

66 

17 

CONTINUE 

67 


01HUX = (3 9HT(X fJ )*U(T ,J,K )«HT(I*2f Jl*UfX-2t JfKI«4*MTfX-lt J)» 

SB 


CU(i-lfJfKn/(2.*OX) 

69 


01HVY = J3*HT(X ,J »*V( I ,J,k»*HT(I,J-2)*v<Ji J-2fK)-4i*HTCX,J-l)* 

70 


CVdf J-lf K))/( 2**DY» 

71 

24 

CONTINUE 

7Z 


RETURN 

73 


END 
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7.1.14 WVEL2 

This subroutine calculates the differential time, in diff- 
erenced form, in the indefinite integral for the equivalent 
vertical velocity at each time step, u, v, and HT are used at 
t - At and D, E, HIE are used at time level n+1, thereafter. 
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C******«***«>4>* *••«****•»* *»•»*•*••»«**»***••«**•*««****•*««*•«*****•••***•>« 

C THIS subroutine CALCULATES THE OXFFErENtIAL (O irFERENCCD) TERHS IN 
C THE INDEFINITE INTEGRAL FOR THE EOUIVaLENT VERTICAL VELOCITY 

C AT each time step 

c*****»«**«*»****«**«*««»*»*»****«******««»*«*****»*»*«*»*«*«***»*»*»*«<»»».» 

SUbROUTINE UVEL2(I,J/ ,IN ,UN ,KN tU tV «HT tOX fOV tHAR ,D1HUX, 0 IHVY, 


COlHUXltOlHVYll 

DIMENSION U(IN* JNtKNliVf IN* J n»KN) «HT( lNt JN) f HAR( IN,JN) 

IF(HAR(I *J|.C0.3I 60 TO 112 
XFIHARiltJI.EO.SI CO TO 119 
IF(MaR(I*U).E0*2) Go To 113 
1F(MAR(I*J) .EO.l) 60 TO 120 
IF(HAR(I«U).E0*4> go to 114 
IFlMARd »J) .E0.7) 60 TO 115 
IF(MaR(1»J)*E0.9I Go To 116 
IFtKARf I*J) .EO.IO) 60 TO 117 

112 CONTINUE 

01HUXr|4uHT(I«l *J )*U Cl «1 * J« K ) - 3 *HT 1 X(J) 4>U (I tUfK) -HT(X« 2* J) * 
CUCI«2*JtK))/{2.*0x) 

OIHUXUI 4*HTCI« I, j)*UCl<»l,J,K-l»-3*HT(l*J)*UClf J,K-1 )-HTcI« 2,J)* 
CUCl«2»U*K.l)|/f 2.*0Xt 

D1Hvy=<HTCI,J* 1 )*V(I ,J4l,K»-HT(l,J-lH>VCI,J-i,K) l/ICt'^OYI 
OlHVYiriHTII,J4 1 )*V C I , J« 1 ,K-1 1 -HT (1 1 J-1 ) *V C I , J-1 »K-1) ) / t 2 .*0 Y I 
60 TO 200 
114 CONTINUE 

01HUX=(3*HT(I *J )*U(I *J *K ) «HT (1-2 *U ( X-2 tJ »K > -4*HT( I -1» J) * 
CU(I-1,J*KM/(2.*0x) 

OIHUXISC 3»HTCI,J)*U(I»J,K-1 l♦HT(I-2 , J>*U ( 1-2 , J,K -1» -4 *HT ( I-l , J > * 
CU(I-l*J,K.in/( 2 .* 0 X) 

01Hvy=<HT(I,J* 1 »♦V(I ,J4l ,K)-HT (I ,J-U*V (I tJ-ltK) »/(2.*0Y) 
01HVY1 S(hTCi,j* 1)*WCI,J*1,K-1)-HT(1 .J'H*V(1*J-1 fK-ln/(2**0V) 

GO TO 200 

113 CONTINUE 

1F(J.C0.1.ANO.I.6C.31.ANO.I.LE.33I GO TO 311 
01HUX=(HT(I4l,JI*U(I'*l»J»K»-HT(I-l , J|*U (l-l ,J,K> l/(2.*0X) 
OIHUXISCHTCI* 1,JI*UII>»1,J,K-1 »-HT ( I -I , J ) *U ( I -I , J ,K - IJ) / ( 2 . *DX ) 
01HVY=(4*HT(I *J41 )*V(I tJ-»l»K)-3<t<HT ( I , J ) 9V ( 1 , J »K ) -HTCltU* 21* 

CVCl* J42* K))/( 2.00Y) 

01HVY1 = (4*HTCI, j4n*V(I,J-»l ,K-1 >-3»HT<ItJI*V(I»JiK-l) -HT(I,J42)* 
CVC:,J42,K-1)I/(2.*0Y> 

GO TO 200 
311 CONTINUE 

01HUX:(HT(l4l,J i»U(l4l,J,K}.HT(I-:*J)*U(X-l*JtK) )/(2 .*dX) 
OIHUXISIHTCI* 1* J»*U(I*1» JfK-1 J-HT (I-l,JJ*u« 2 .* 0 X) 
01HVY=(4*HT(I ,J4U«V(I ,j4l,t<)-3«HT ( I , J I * y ( I » J »K » -HT(I» J« 2 J* 
CV(I,J42,K))/( 2.«0Y> 

OIHVYUC 4*HTCIt J*l)*V(I,J4l ,K-n-3*HT(I,JI«V(I,J,K-l) -HTC I, J42)* 
CV(I»J424K-in/(2,*0Y> 

60 TO 200 
120 CONTINUE 

IF|J.E0.11.AN0.I.6E.7.AN0,I,LE,16) GOTO 321 
DlhUXr|HTCI*l,J HiUCIM ,JfK)-HT(I-l , J)*J (I-l , J ,K ) ) / ( 2, *0X I 

01HUX1=IHT(14 1, jj *UCI«1, j,K-l >-HT (l-l , J)*u< I-1»J *«-!) )/ C 2.*0X ) 

01HVY=(3*HT(I ,J J*V(I ,J,K >4HT(I,J-2) »VCI,J-2,M-4*hTC 1, J-l»* 

CVCIt J-l*K»>/( 2.<'DY) 

D1hVY 1=C 3*4T( I, J)t'V(I*.J,K-l »-*HT» ItJ-2)*V a, J-2,K-1>-4«HT{I,J-1I* 
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CV(ItJ-lf K>l>k/<2.»0Y) 

CO TO 200 
321 CONTINUE 

01HUXr(HT(l«l ,J l»U(I«l (J ) .HT(I>t» J7«U(I*t *J»KI l/(2.*0X> 
OlHUXlstHT(Z« It J)»U«l«ltUtK-l )-HT (1 >1 » J t *U ( I •! » J t« -i»l / ( 2 **0X I 
01HVVr<3*HTCI ,J i*V(I tJiK I ♦HT ( I , J-2 > *V < 1 1 J-2 »K I -4 *HT (I t J-1 > • 
CVfltU~ltKn/( 2.<»0YI 

OIHVYUI 3*HT(lt%n*V|I ,J,K-1 I ♦HT ( I ,J -2 1 *Y < 1 1 J-2tK -U -4«HT i 1 1 J-1 » • 
CV41tJ'lrK>*l. )/;2t*DY) 

CO TO 2GC 
ns CONTINUE 

01HUXr(4«HT(I<»l )»U (I «1 , J,K )- 3 >«HT ( Z tJI*U 1 1 , J»KI -HT(I«2t J7* 
CU(I«2f JtKM/(2*«OX> 

OlHUXm 4#HT4 !♦ ltUl«UCl«l tJtK'l )-34>HT(I tU)*U(Zt JtX>l|*HT(I«2, J)* 
CUII«2tJtK>in/(2.*0X| 

01HVY=(4<»HT(I tU«l )«V (I ,J«1,K )>3*HT ( 1 1 J I *V ( I t J tK I 'HTl I « J« 2 ) * 
CV(ItJ«2tKn/(2.*0Y) 

OlHVYl=l 4*HT(I, J^lI^vatJ-*! fK-1 l-3*HT<It j7*VCIf JfK-l» -HT<It J«2»* 
CVdf J«2tK>l||/(2.*0V| 

GO To 200 
119 CONTINUE 

DlHUX=444HT«I*i tJT^UJI^ltJfKJ-J^HTi I,J7*U(I,JtK» -HT(I*2t Jl* 
Cu(l«2tJtKn/(2.*0X) 

01h>^X1 = ( 4*hTI ltJI«U(I«ltUtK>17-3>»HT(ItJ »*U<I,J,K-1» -HT(I«2,Jl* 

CU( I«2,Jf K-1) ) /( 2.*0X) 

0IHVY=(3*HT(I «J tJtK » ♦HT C 1 1 J-2 I *v< I » 0*2 t K I -4 *HT4I t J- 1 )* 

CVdtJ-lf KII/4 2.*0YI 

01HVY1 = < 3*HT( I, J)*V 41 ,JtK-l WMT{I,J-2l*va,J-2tK -1 1 -4*HT 4 1 1 J-H • 
CV41,J-1 ,k-1I»/42.*OYI 
60 TO 200 

116 CONTINUE 

01HUX = C3*HT4I ,J l*U <1 ,J ,K I *HT 4 1 -2 , Jl '“U 4 1 -2 , J ,K ) -4*HT 4 1 -1 1 J > * 
CU4l>l«JtK)t/42.*OX) 

01HUX1 = 4 3*HT4 It J»*U<I,J|K-1 l«HT4 1-2 , J l*U 4 1-2 1 Jt*< -1 3 -4*HT ( I-l, Jl* 
CU4l-ltd|K-in /4 2,*OXI 

OlHVT=44*HTCI ,J jl tJ«l t*< I'^^HTi I,J1*V 41 , JfK) -HT4I tU»2>* 

CV4ItU«2tKM/42.*OYI 

01hVYl=4 4*HT4It J^IH»V4I t J-»l tK-U-3*HT4ItU»«‘V( It JtK-ll -HTtIt J«2I* 
CV4ItJ*2tK-l)»/4 2.*0Y| 

60 TO 200 

117 CONTINUE 

01HUX = 4 3*HT4I fJ l*U4I ,J tK » •*H T 4 1-2 t *U 4 1-2 t J tK I -4 *hT< I-l t J> * 
CU41-lt JtK)»/( 2. 90X) 

DIhUX1 = 4 3*hT4 1 1 J» *U 4 1 , J tK -I ) *HT 4 1 -2 , J) *U 4 1-2 1 J t>< *1 » "'♦♦HT 4 1-lfU)* 
CU4I-ltJtK-l)l/42.*OX» 

DlHVVr(3'»HT4I ,J »* V4 1 , J ,K ) ♦HT 4 1 1 J- 2 1 *V 4 1 1 J -2 t K » -4 i»HT 41 t«i-H* 
CV41,J-ltKH/42.*OY» 

01MVY1 = 4 3*HT4 It J30V 41 ,JrK-lMHT 41iJ-2»*V 4 1, J-2,K -1 J -4*HT4 I, J-n* 
CV4I,J-ltK-l) ) /4 2«*DY) 

200 continue 

RETURN 

ENU 
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7.1.15 WW 

This subroutine converts ^2 to W, that is the vertical 
component of velocity. The following analytic expression Is used 
for this conversion: 

W - AH + (c9^ + (0-1) ^ 

The actual vertical velocity component, W Is defined as WZ(X,J,K) 
In the model program, and It Is calculated at each x, y,o. 

Since WZ(I,J,K) Is not used In solving the system of governing 
equations, this subroutine Is used only after the last time 
cycle for each computer run. 
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THIS SUBROUTINE TRANSFORMS THE EQUIVALENT VERTICAL VELOCITY 
(IN THE SIOMA COORDINATE SVSTEMI INTO THE ACTUAL VERTICAL VELOCITY 
(IN the X-Y-2 coordinate SYSTEH| at each X-Y LOCATION AND DEPTH 
IN THE domain 


SUBROUTINE WU (I Nt JN ,KN ,H T D , HT E * ET A , D tE t V »UZ »M AR , DX «CY , DZ t DT t HX , H Y ) 
dimension HTD(lNtJN)|HTE (IN .JNI .ETA ( IN , JN I . W ( IN , JN .KN ) . 

CUZdNf JN«KN).MA:!(IN.JN) tD(IN.JN.KN) .EdN, JN.KN) tHX(IN.JN) 
C.MYdN.JNI 
KNl=KN-l 
DO 10 1=1. IN 
DO 10 Jrl.JN 
DO 9 K=1.KN1 

QZ(ltJiK>=0»0 

ir(MAR(I*JI.EQ«C> GO To 9 
IF(MAR(I .Jl .EQ. 11) GO TO 11 
IFlHARd.JI.EC.C) GO TO 11 
IF(KAR( 1. J) .EO.F) GO TO 11 
IF(MARd.J) .EC. II CO To Icl 
IF(MARd.J) .EQ.?I 60 TO 102 
IF(HAR(1 .JI.EC.Z) GO TO 1C3 
IF(MARd.J).EC.u» SO TO 1C4 
IFlHARd.JI .EO.S) GO TO 105 
1F(MAR(I.JI.E0.7I GO TO 1P7 
IF(MAR(I.JI.E0.9) 60 TO 109 
IF(MARd.JI.EO.lO) GO TO 110 
CONTINUE 

ETAX = (ETAII4 1.J)-ETA(I-1,J) )/{2,*oX » 

ETAY = |ETA|I,J*11-ETA(I .J.IM / (2.9DY I 
60 To IflO 

ETAX = (£TA(I*1.J )-€TA(I-l. j) )/(2.*0X ) 

ETAY=(3*ETA(1 .J )*ETA(I . J-2 » -4»ET A ( I .J-1 ) I / ( 2 .*0Y I 
CO To 100 

ETAX = (ETA(I*l.J l-ETA(I-l.J) )/(2.*0X > 

ETAY=(4*ETAd .J ♦ 1 )-3*E T A ( 1 . J » -ET A ( 1 . J« 2 » ) / ( 2 .*0Y ) 

GO TO 100 

ETAx = (4*ETAd4l ,J}-3*ETA (I, J)-ETA(I ♦2.J))/( 2 .*DX ) 

ETAY = (£TA(I.J«1 l-ETA(I.J-l) l/(2»*DY ) 

GO TO ICO 

LTAX = (3*ETA(I .J I ♦ ET A (I -2 , J) -4*E T A (I -1 . J )) / ( 2 .*0X I 
ETAY = (£TA(I.J^l)-ETAa ,J-1U/(2.*0Y ) 

CO TO 100 

ETAX = (4*ETA(I41 ,JI-3*ETA d, J)-ETACI ♦2.J))/( 2.*oX) 

ETAY = (3-»ETAd .J )*ETA (I ,J-rj-4*ETA (I .J-n ) / ( 2.*DY I 
GO TO 100 

LT;.X = (4*ETA* I ♦! .J)-7*ETa 1 1. J)-CT A( 1*2. J) ) / (2.*0X ) 

ETAY = (4*ETA(I .J^n-I^ETA ( I, J)-ET A(I ,J*2) ) / ( 2.*0Y » 

Go To iQC 

ETAX = (3*ETa(I ,J »*ETA(I-:,J)-4*ETA(I -l.J) I/(2.*0XI 
ET/.Y = (4*ETA(I .J«1)-3*ETA (I, JJ-ET a(I ,J«2))/12.*0Y I 
GO 10 100 

ETAx = (3<tETA( I .J ).ETA |I -2 , J) -4*ET A(I -1, Jl I / ( 2.*DX » 

ETAY:(3*ETA( I .J METACI , J-2)-4»ET A (I ,J-1 l|/ (2 .*DY I 
COM INUE 

(I. J,K )=HTE (I ,J )*W(I ,J ,K ) ♦ ( (K -1 )*( 01TE( I ,J)-HTo( I . JU / PT 
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57 C«OU |JtK METAX«r(I»J»K I*ETAVI«( IK-1 l»DZ )• CD I ItJt Kl 

58 C*HX|ItJMCClf<*'<|*HVlX,J>) 

59 9 continue 

83 ID CONTINUE 

81 RETURN 

82 CKO 


T "CTFntM'^ pi 1,^.,, rnr nj jupp) 
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7.1.16 PRES 

This subroutine calculatts the pressure field at time level 
n+1 by knowing the contour dppth, H and density, p at n+1. Note, 
that this is the integrated form of the hydrostatic equation. 

The integration is performed by applying the trapezoidal rule. 
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1 c 

2 C************* ***«*••**•**«***•******«*••• -«**•******•*•«***•«>•*•**» 

S C this SunKCUTXNC CaCULATCS THC PRCSSUOC FXCLO 

% c************* ****•**«'•**«**«•**•'»•»•«•**'»*'•»»•*••*•*•**»*****••*•*•• 

5 SUURCUTZNC PRF.S IlN, jN»kN ,HT »RO tCH tP tCZ > 

I OZHFNSION NTIlNtJN>tRO(ZN»JNt><Nt(PIIN»JN(KNI 

7 00 1C I=t*ZN 

• 00 10 J=1(JN 

9 P(Z(J»1I=3.3 

19 00 e K=2tKN 

11 P(Z,J,K)=P(I,J,K<-1MCR*HT(X »J1*IR0( ZtJ(K-l)«ROf X iJfKI I *92/2 .0 

12 • CONTXHUC 

13 10 CONTINUE 

19 RETURN 

IS tNO 


"'’''A'. r- 

#'UOH quality 
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7.1.17 TEMP 

This subrouclne calculacer Che cemperacure discrlbuclon 
T(I,J,K) ac X, y,a at each cimc step using a forward differencing 
in cime. T(x, y, a) at C- 4t (■T(1,J,K) is calculated from 
T(x,y,o) at t"0 (-TNICI.J.K)) . Thereafter, T at time level 
n+1 (-TN1(I,J,K)) is calculated from T at time level n(-T(I , J ,K) ) . 
The spatial derivatives, once again, have been approAl....'.*-'''^ by 
central differencing in the interior of the domain, and three 
point single sided differencing on the boundaries, except at 
MAR"6 and MAR«8. The numerical scheme is given in Volume I. 

Note, that Che adiabatic approximatio i given in Volume I 
calculates Che temperature at boundary points after tne energy 
equation is solved at interior points and at MAR«6 and MAR*8 . 
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23 
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29 
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33 

31 

32 

33 

34 

35 

36 

37 

30 
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C*****»*««***«»*4i***«»***»»**»«*«*«*****«**M9**«**«9«*»*»**99**4i***<»*»«*« 

C THIS SUBROUTINE CALCULATES THE TEMPERATURE DISTRIBUTION AT EACH X-V 

C LOCATION and DEPTH IN THE DOMAIN AT EACH TIME STEP USIN5 A 
C rOfiWARO OirFERENcING SCHEME IN TIME 

C»»*****W«»«««»*«**«»»*»9»»»*****«i*«»*«********«**»«***«*«*««****<l>V****»** 

SUBPOUTINF TEMP (IN,JN,KN ,T,U,V,W,DX »Cr .OZ tOT ,H,3H»DV, MAR, TN I,HN1 » 


CORAO,RR,HS, TA.ROI 

DIMENSION T(ZN,JN,KN),U(IN,JN,KN>,V ( IN* JN ,K N > , W ( IN , JN ,KN) ,H ( I N, JN 1 
C»MAR(IN, JN),TN1 (IN , JN ,KN ) ,HN 1 ( IN * JN ) 

C*ROf IN,JN,KN) 

DO 10 I:1*1N 
DO 10 USltJN 

IFtMAR(I,J|.EO.Ot GO TO 10 
DO 8 KS2,KN 

IFIHARI I ,JI «EC. 11) 60 TO 11 
1F|MAR|I,J).E0.6) SO TO 11 
IF(MARII,J).EC«e) SO TO 11 
IF(HaRCI,J) .LT. 11 ) GO To lO 
11 CONTINUE 

OHXr(H(I*l,J)-H a>l ,J) )/ (2* mOX) 

0HY=(H( I*J«1) >H (I ,J-1) ) / IZ.O'OY) 

01TX = m I<*l,d,K )-T(I-l ,J ,K) )/(2.*0X 1 
01TY=(T( I,J«1,K >.r(ItJ-l,K) )/( 2 *« 0 Y) 

02TX=CT( I-»1,J,K )^T(I-1 ,J,K)-2.*T|I,J,K))/{OX*OX> 

OZTYsm I,J*1,K •♦TJI,J-1 ,K)-2.*T(I. J,K) ) / (DY*OY ) 

01HUTX=IH(I«1,J )»U(I*1,J,K)-«T(I*1,J,K)-H(I‘*1 1 J)*U(I-lt J»K)* 

CTf I>l»J*K))/f 2. «0X) 

D1HVTY = (HII,J«1 )x<V(I,J«l tKIOTdfJ^l ,K)-H(I,J>1)*V(I,J-1,K)9 
CTa,J>l,K))/(2.*0Y) 

IFiK.EU.S) 60 TO 71 

01HT2=(W«I,J,K4 1)<.T(I,J,K*1 i-W(I,J,K- 1 )*TtI,J,K-l))/( 2,*02I 
02TZz(T( 1,J,K*1 )<TII,j,k-1)-2.*T(I, J,K))/C0Z*02) 

60 TO 80 

71 0lHTZr|3*TCl , J, K » ( I ,d ,K )4T(I (I ,d ,K-2) -4*T( I,d,K-H* 

CN(I,J,K-1))/(2*M0Z) 

02TZ=2.*(TI1,J,K-1)-T(1,J,K »)/ (D2*0Z> 
aO CONTINUE 

TC=<D1HUTX*D1HVTY4H d,J)*Dl HTZ) 

TKX=BH*I0HX*0 1TX4H(I ,J J*D2TX) 

TKYrBH*(OHY*0 IT Y*H< I,J)*02TY» 

TK2rBVM02T2/H( I,d) ) 

TN1(I,J,K)Z(I -TC + TKX4TKY4TKZM'0T4H< I,J)*T(I ,J,K) )/HNl(I,d» 

8 CONTINUE 

10 CONTINUE 

00 1000 IrldN 
00 1000 JZl,JM 
IF«MARIi,j).E0.5J go To 1030 
IFCHARIItJ) ,EQ,6) GO TO ICOO 
IF(MAR«I ,J) .EC. 7) 60 TO 1000 
IF(HAr(I,JI.EO.Q> go to loco 

IF(MAR< I f J) .EO. C) GO To ICOQ 
IFCMARI I,J) .EC. 10 ) GO TO ICOO 
IF(MAR| I ,J) .EO. 11 ) GO TO ICCO 
00 lC3l K=2,KN 

Ir(MAR(I,J) .EO.C) TNl (I ,J,K )=a.O 
imurni.j) .EO.l ) TNICI ,J,K >=TN1 (I, J-1,K) 
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57 


lF|HA(7t It TNlIXtJfK >:TN1 ( t, J«1 1 K 1 

SB 


IF«KAR(I,J}.eC.3l TNllItJtK)STNl (X«ltJtKl 

59 


Xrcr'ARd tJ»tEC.9l TN i(ItJ»K)=TNl<I'lt«7»Kl 

50 

1001 

CONtI^^UE 

51 

1000 

CONTINUE 

52 


00 2003 XsltlN 

53 

' 

00 2000 JdtUN 

59 


XptMARC Itjl tCC.C) CO TO 2000 

55 


IFtHARd tJI.EQ. n CO TO 2000 

56 


XF(MAR( Xt J> *e0.2l CO TO 2000 

57 


XFtMARd tJdCG. 3) 60 TO 2000 

5B 


XF(MARdtJ)*EQ.'il SO TO 2000 

59 


XFIMARdtJl tCO.O CO TO 2CS0 

73 


IF(MA«d iJ» .EC.e 1 Go TO 20C0 

71 


XFtMARCItJl.EO. Ill 60 TO 2000 

72 


DO 2001 K:2,KN 

73 


XF(NARd tJI.EC.S) TNKIyJ.K ) = (TNllX«ltJtK>«TNldtJ-I|K} 1/2. 

79 


iFtHARt It JI.EO. 7) TNllltUtK I = (TN1 II41,J,K H TM (I yj4l,K}|/2. 

75 


IRKARd tJI .EO.9) TNKI t JtK l = (TNllI-ltJtK >«Tk 1(I t J-*lt KII/2. 

75 


lF<NA9d tJl.CO. lu I TNKI tJt*< > = (TNi< I-l, J t K I 4 TN K 1 1 J-1 t K I 1 /2 

77 

2C01 

CONTImUE 

78 

2000 

CONTINUE 

79 


DO 100 I=ltIN 

BO 


00 100 J=ltJN 

el 


If(MAR( ItJl tEQ«CI TNKItJtl 1=0*0 

B2 


ROdfit 1 1 =1*029 93I-.COC02C*TdtJtl)'*0000098*(Ti ‘,1I«*2I 

B3 


RRzROfX, Jyl) 

B9 


TERN=(D2<»MSitHNl II ,J )/ (RR4>BV II 

65 


TNldtJt ll=CTNl(I,Jt2l*TA*Tt:RN>/(1.04TERM) 

B6 

100 

CONTINUE 

B7 


RETURN 

BB 


END 


PR£C£DIM'^ 
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7.1.18 OLDHT 

This subroucine transforms H(x,y,o) at time level n+1 
(-HTE(I , J)) into H at time level n (-HTD(I , J)} and H at time 
level n is transformed to H at time level n-l(«HT(I , J) ) . These 
transformations are performed at the end of each time step. 
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7.1.19 OLDUV 

This subroutine tranforms u and v at time level n«l 
(«Dv.I,J,K) and E(I,J,K)) into u and v at time level n(»H(I,J,K) 
and G(I,J,K)), and u and v at time level n is transformed into 
u and V at time level n-1 (“U(I,J,K) and V(I,J,K)). These trans- 
formations are performed at the end of each time cycle. 
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1 

2 

C 



3 

C 

THIS SUBROUTIK'E TRANSFORMS MATRICES 0|E INTO U,V| BCSPECTIVELV, 



c 

FOR the next time CtCLE CALCULATION OF HORIZONTAL VELOCITIES AT 


S 

c 

EACH X-V LOCATION aNO CEPTh IN ThE DOMAIN 


6 


7 


SUUROUTINE OLDUViIN ,JN ,KN ,U vVtHtC»0 tCl 


8 


DIMENSION U(1N,UN»KN |,V(1N,JN»KN) •D(INtJNiKN) tE I 1N,JN,KNI t 


9 


CH( 1N,JN,KN> »8<lN»UNtKN) 


10 


DO 10 KS1,KN 


11 


00 10 iSliIN 


12 


00 10 J=1»JN 


13 


U(l,J»K)rH(I,J(K) 


r<i 


V(I»J»K)=G(I»J»IO 


15 


H|1,J,K)=DCI*J»K) 


16 


C(I, J,K)rE(I,J,K) 


17 

10 

CONTINUE 


18 


RETURN 


19 


END 
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7.1.20 PLOT 

This subroutine transforms T at time level n+1 (-TN1(I ,J ,K) ) 
into T at time level n(-T(I , J,K)). This transformation is per- 
formed at the end of each time cycle. 


o n o o 
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THIS SU8R0UT1KL TRANSFOPHS HaTOIX TM INTO T FOR THE NfXT TIME CyCUE 
CALCULATION OF TEMPERATURE AT EACH X-V LOCATION AND 3EPTH IN THr oOH/. 

SUOROUTINE OLOT (INtUNtKN •T»TN1I 
DIMENSION T(lN,JN»KN|,rNl(IN,JN,KN) 

DO 10 K:1,KN 
00 10 isltlN 
DO 10 JsltJN 
TII«J»K»:TNm» J,K1 
10 CONTINUE 
RETURN 
END 
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7.1.21 ETT 

This subroutine calculates the wave neight 
(■E1A(I,J)) at the end of each time step. 

ETA(I.J) - HT(I.J) - HI(I,J) 


t 


n (x,y) 


o o o n o 


THIS SUSnOUTXKC CALCULATES THE WAVE HCIGHTISURrACE ELEVATION ABOVE 
HEAN water LEVEL) AT EACH X-Y LOCATION IN THE DOHAIN AT EACH T1*1C STL 

**•***«•<»••«• ««>««*** •****««***«****««*••**¥*«*** ****A*^t ««**•*»«*»«** 

SUBROUTINE E T T t IN »JN tHT tHI |H AR tCT A ) 

DIMENSION HT(INfJN)(HI (IN,JN),NAR (1 N ,JN ) , ET A ( IN* JN) 

DO 10 I=ltIN 
DO 10 J=ltJN 

ir(HAR(I»J|.EO.C) CO TO 1C 
ETAt It J)=HT(IfJ)-HI(I(J) 

10 CONTINUE 
RETURN 
ENO 
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7.1.22 PRPARA 

This subroutine writes out the physical and numerical 
parameters at the end o£ each computer run, i.e. after the last 
time cycle for a particular computer run. 
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1 c 

2 C»**»*»**»***»0**»**0f*0*0*»***0*0»00m4^*00*0»0t»*0*t*»**t0t»»»0*»0*»04tt 

3 C IMIS SUBROUTINE PRIMS TMf PHYSICAL AND NUMERICAL PAB*METERS rOR Int 

1 C variable PENMTY model at TME ENU or EACH COMPUTER RUN 

5 C *****•••«*•«• 

fc SUttROUTINf PRPARAICI,CH,CV,CP,CC,OX »OY ,PE ,0 T, T*UX ,TaU Y, TTOT ,6Rtrr • 

7 CRR,NHtKv,bH,0V»ORAO»TI ,T TOT 1 ,T A ,HS 1 

• IF«TTOU.cT.O*l ITOTIsTTon-OT 

♦ PRINT ltCI,CH,CV,CC,CP,OX,OY,Or,OT,TAUX,T AUY,TTOT,CP»rr,HR,MH,Kw 

3 C»eH,PV»ORAO»TI(TTOTI,rA,NS 

1 1 FOKMATC/* CX = *,El!i.7,/» CHs • ,C15 . 7 , /• CVs*,CI5.7/« CP = »,E15.7, 

2 C/» CCS*,EI5.7»/* 0Xr*,L15.7,/» OVS* ,E15 .7 ,/ • OZr * ,EIS . 7 , / * OTr», 

3 CE15.7,/* TAUXS* ,El5.7,/» TAUYs* ,EI$ .7t/* TTOTr*,E15.7f/* gR=*» 

V CEI5.7*/* FFS» ,E 15.7,/* RR s* ,E15 . 7 ,/ • KMr«, EI5.7,/* KWr*,Ei5,7/, 

5 C/* QHr%CI5.7«/* BVr * ,E 15.7 ,/ • OR AOsSElS .7 ,/ • TX=*,£15.7|/ 

6 C* T70Tls»,E15.7,/» TA = », £15 . 7 t/ * MSs*,E15.7/l 

7 ttoti=ttoti*ot 

9 RETURN 

9 END 
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7.1.23 PRETA 

This tubroutlns wrlcct out ETA(1,J) «t the end of the computer 


ruQ. 



1 c 

2 C************************************************************************* 

S C THIS SUBRCUTXKE PTINTS THE WAVE HEIGHT AT EACH X -V L3CAT10N IN THE 

% c oohain at the eno or each computer run 

S SUBROUTINE PRCT A ( I , J ,IN , JN« ET A ) 

S DIMENSION ETAIIN»JNI 

T DO 10 I:l,XN 

S 10 PRINT 11 (ItfCTA (I ,J| yjzl ,JN ) 

V II FORMAT!/* I :*|1S/* UA VE-HC ICHT '/ ISXtAElS .7 1 1 

19 return 

11 END 
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7.. 1.24 PRUV 

This subroutine writes out U(I.J.K) and V(I,J,K) at the end 
of the computer run. 


o o o o 


»*«***»* 4^**»*W*<*»4<**|^«***»***^*««» «*•#*•** **»»•***•» *•«***•***»***« 

This suu(<outinc points thc horizontal vclocitils u,v at lacm x-v loc 

AND DEPTH IN THE pOHAiN AT THE EnO OF EACH COMPUTER RUN 

*******<^*0>!i**#i>*i^i|i*4»4M>***0 0*0***W4>«<|i4i«***<|i>t>****4>0***#«^0*4i4>0#$4i**«*>' 

SUBROUTINE PRUV (I *J ,K JN , JN ,KN ,U • V ) 
dimension UIIN» JN*KN>*V( IN» JNfKN) 

KN1:KN-1 
00 10 K:1,KNI 
00 10 I=ltlN 

PRINT 11 tK«I» (U (1 »J«KI »jrl, jN) 

10 PRINT i:,(Vn«JtKI»J = l»JN> 

11 FORMAT!/* Kr*,I3,3X,«l:*,I3,/* U -VE LOCITT ♦/ ! sX 1 6 E 15 .7 ) ) 

12 rOKMATI* V-VELOClTY*/l5X,0£l5.7n 
RETURN 

END 
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7.1.25 PRW 

This subroutine writes out WZ(I,J,K) at the end of the computer 


run. 


1 C****************»****************>^***********'0*****m********************«>f 

Z C THIS SU3R0UTINE PRINTS THE ACTUAL VERTICAL VlLOCjTY AT EACH X -V 

5 C LOCATION AND DEPTH IN THE DOMAIN AT THE END OF EACH COMPUTER RUN 

A C** *********"'******************** *"****«4'*'AiAW»**<AA 4*** ***<A*A*****'A*«>4i» 

5 SUBROUTINE PR W( IN »JN «KN »!• ) 

6 DIMENSION WIINtJNfKNI 

r KN1=MN-I 

S DO 10 ’^SlfKNI 

9 DO 10 i:l»lN 

10 10 PRINT ll»K*I»iU jNI 

11 11 FORMAT!/ • K = » , I J,3X ,* 1 =*,l3|/» U -VELOCI TY •/ ( SX , BE 15 . 7 » I 

12 return 

13 END 
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7.1.25 PRTEMP 

This subroutine writes out T(I,J,K) at the end of the computer 


run. 


1 

2 

3 

d 

5 

6 
T 
B 
9 

13 

11 

12 

13 


C THIS SUBROUTINE PRINTS THE TEMPERATURE AT EACH X -Y LOCATION AND DEPTh 

C IN the domain AT THE END OF EACH COMPUTER RUN 


SUbROUTINE PRTEMPd ,J,K tIN| JNvKNiTI 
OIHENSION KINfJNtKNI 
DO 10 K=1«KN 
DO 10 I=ltIN 

PRINT lltKtIf IT(XtJ,KltJ=lt JNI 

11 FORMAT!/* K = *,13,3Xt* I=*tI3,/* T EMPER ATURE * / ( 5X t^E 15 •? I ) 
10 CONTINUE 
RETURN 
END 
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7.1.27 STORE 

This subroutine writes on magnetic tape all calculated 
system variables and numerical parameters, DX, DY, DZ, DT, TTOT 
and TTOTl. 
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1 c 

2 C************* »••♦♦*•■»■»*♦•*•••••••■••••»*• *•*♦••••••••••••••••♦•**•*••*•*• 

3 C THIS subroutine; writes on hacnctic tape* for the variable density 

R c houel* the values for the variables and physical and numerical 

5 c parameters for storage and for reading in data for the next compute 

6 

7 SUDROUTINE STORE (IN ,JN ,KN ,U ,V ,M ,HI ,HT tHTD »HX ,HY ,HAR,ETA,P,RO,Clt 

I CCCiCMtCV,CPtOX,rY.OrtOT.TAuXtTAUYtTTOT,H,6,HTE,T,TTOTl,M2l 

9 DIMENSION U(I Nt «^N »KN I ,V (IN , JN »KN 1 ,U ( jN » jN ,knI »P< IN>JN(KN) » 

13 CHIIIN«JN|.HT(IN,JN) ,HTD(IN, JN},HX(1N,JNI,HY lINtJNI tH«R(IN,JNIt 

11 CETAf lNtJNltRO(lNtUN(KN)»H(INf JN»MN) »C( IN , JN ,HNI , HTE ( IN, JNI 

12 C,TfIN,JN,KNI,W2(XN,JN,KNI 

IS WRITE (6MIIUII,J,Kl,Ksi,KN),J=l,JN),lri,INI, 

IV Cmv(l,J,K),Krl,KN),j7l,JN) ,lsl,INI , 

15 CKtWIi, j,K)tK=l,KN) ,U:1, JNI , i:l ,1N) , 

16 CC((H(I,J,KI,K:l,KN),J:l,JN),x = l,lNI , 

17 C( ((6(I,J,K),KSl,KN) ,J=1, JN) ,l;i,IN) , 

19 C(Clp( I, J,K},Kri,KN) fjrl, JNI tIslyIN) , 

19 CdIROd, J,KI,K=I,KNI,j7l, JN t,l = l,IN>, 

2 3 C( iHTDd, JN),l7l,INi, 

21 C< IHlEd, JN),I=1,IN I, 

22 CUHId,JI,J71,jN),I71,lNl« 

23 Cl (HXd, J) ,J=lfJM ,l=l,INt , 

2V C((HY(I,J),J=l,JM,l7xtIN)f 

25 CHMARd, JN),I=1,INI, 

26 C((HTd,JI,J=l,JNI,Irx,lNI, 

27 CUETAd, JI,J=l,JNi,I=l,IN|, 

29 CM(Td,J,K>,Krdt^NI,J=l, UN) tl^ltlNI , 

2 9 CdlUZd, J,KI,K=1,KNI ,J = 1,JN l,I71,INI, 

S3 CCI,CC,CH,CV,CP, OX,Oy,OZ,OT,TAUX,TAUY,TTOT ,TTOTl 

31 END FILE 8 

32 RETURN 

S3 END 
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7.2 MAIN PROGRAM FOR NASUM III 
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7.2.1. TMAIN3 

This is the main program for free-surface complete 
field model. This program reads in the data, initializes 
the necessary quantities, co-ordinates the subroutines and 
calculates the velocity and temperatures in the whole domain 
under consideration. The parameter statement defines the 
size of the computational domain. The subroutine "XYSH" 
does horizontal stretching. The subroutine "READ 2" reads 
the MAR matrlcs which distinguishes the various points in 
the domain. The subroutine "INITIB" sets the initial condi- 
tions on velocities, temperatures, surface height and reads 
the depths at various points in the domain. The subroutine 
"CURNT" which is called after "INITIB" sets the velocities 
everywhere in the domain equal to the current velocity. 

The subroutine "INLET" puts the discharge velocity and tem- 
perature at the discharge location. Then it follows a set 
of subroutines to calculate the velocity and temperature 
field for the entire domain. The values of variable at 
different time levels are given in the Table (1) . 
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I • 

1 



I 

? 

3 

4 

5 
i 

7 

8 
9 

10 

n 

12 

13 

15 

16 

17 

18 

19 

20 

21 

22 

21 

2*1 

-»c 
** « 

27 

2 ? 

29 

31’ 

31 

32 
37 
3't 
3«^ 
3^ 

37 

38 

39 

40 

41 

42 
4 3 

4 4 
48 

46 

47 

48 

49 

50 

51 

52 

5 3 
54 
54 
56 


PARAHfTCR INr,20»JNS20|HN=S 
REAL KH»KV|AH,0V 

DIMENSION U(IN,JN,KNI tV( IN,jN,NNI,NIINtJN,KN| ,UM(IN,JN,MNlt 
CHT(TN,JN 1 ,HI ( 1N,JNI .ET A( IN,JN IfMTniN ,JN) «VMI iNt JN.KNI* 

CXXIINi »VY IJNl ,X»X< INI , VYr ( JNI «X IINI,V UN) ,A (INI,»( JNI , 

CH(IN,J^ft'^> *C ( IN iJN,KN 1 »0 UN *JN ,KN) ,£ ( I N , JN *KN I , WH UN « JN , KNI , 
CMARIIN, JNI,H1DUN, JNI tPO(lN,JNtXNItP(IN«JN,NNI 
C,T(IN,JN,KNI ,TN(IN,JN,KN ) tlF UN»JN,KN),TAMIIN,JN(NNI 
READ ItIRUN 
READ I,LN 
READ i.KSTcrrf 

1 format 1151 

READ 2« Cl,CCfCP.CH,CV 

READ 2| GPtFF,RR,MK 

RFAO 2* OX,0Y,02 

READ 2 , KH,KV,BH,BV 

REAP 2,ntLX,0ELY,CtEX,0EEY,ecrX,EEEV 

CALL XYSHUN»JN,OELX,OELY tDCEXiOECYtEECX tCCCYtXX »YY,XXX ,Y VY , A *0* 

CX ,V I 

2 FOW“Ari 1 

TFC IRUN.GT.O I CO TO 4 
CALL >?EAO?( IN, JN »MAPl 

CALL INlTIn f IN.JN «KN,CTA ,H T (HI ,U t V «P0 ,P «GR »PR »0? •tiTO ,0 ,E «H ,G f 
CHTE ,T,TN , Tf , TAMU^tYMtUtOXtOY 1 
CALL CUC'NTl I »J,K , IN»JN fKN ,U|V.H,G«OtC > 

CALL INLETd ,J«KtlN«>iN,KN,WM,r,TN,TF) 

CALL WMTOVM fN«JN ,KNtHTC«HTO*HTOtCTA,HtC»UtttHyOX«OY(OZ,OTt 
CHAR, XX, YY I 
TTOTSO.O 
GO TO 6 

4 CONTINur 

CALL n'EAPK TN, JN,KN,U, V,N,HI,Hr,HTO*HAR,ETA,P,RO,CI,UH,VM, 
CCC,CH,CV,CP,OX ,PY ,02,0T,TAUX ,TAUV,TT0T,M,6,MTE,T ,TN,TF,TAM,TA!R,D, 
CE) 

6 CONTINUE 

RFAP 2,TAUX,TAUY 
REAP 2,TA1R 
READ ’,DT 
00 4 L = 1 ,LN 
TTOT=TTOT«OT 

CALL UVELJUN, JN,KN,U,V,H,C,0,E ,0X ,0 Y ,0 2 , U , T AUX , T AUY , OT , 
fHTO ,HTO ,HTE ,HX ,HY ,ETA,P,MAP,KM,KV,f;R ,RR ,FF, 

CCP,CC,CI ,CH,CV,RO, IN,XX,YV,XXX,YYYI 
CALL 4 VEL (IN ,JN,KN ,H,0 ,W ,HTO,nX ,0Y,P2 ,:<AR ,XX ,YY,WH1 
CALL TENS ( IN, JN,KM,HI0,H TP,UTE ,OX,OY ,CZ,0T,5H,BV ,T ,TN,TF, 
CW,H,G,MAR,HK,TAIR,1AM,K0,XX,YY,XXX,YYY,L,LN1 
CALL 1 0UNC2 ( I, J,K , IN, JM,KN,U,V,H,G,0 ,E ,M,NI ,HTE, T ,TN ,Tr I 
CALL PFNS1 Y ( IN ,JN ,KN,RO, TF 1 
CALL PP^S (IN,JN,KN,HTE ,R0,GR,P,02I 
CALL f. TT (IN, JN,H TF ,HI ,MAP ,LTA 1 
CALL OLPHf(rN,UN,Mtf ,HTP,HTI 
CALL rLPUVT ( IN ,JN,KN,U ,V ,U,G ,0,r ,T,TN ,TF1 
CALL INLETd ,J,i< , 1 4 , JN , K N , WM , T , TN , TF ) 

CALL Ml y T (I , J,K, IfJ ,Jfl ,KN , TNI 

CALI lNLrT(I,J,K,Ii*,JNfKN,4«l,T,TN,TP 1 

CALI. T IPE (1 , J,K , IN ,JI, ,«N ,ll,v ,M,f. ,n,F , r,TN ,in 


i 


/ 
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•,7 'j CONTIKUr 

sfi iriKSToor.cT.r ICO to icoo 

i.<t CALI. ! TOKf (TN»JStNI«»U»V»W«Nl»HT »HTO t »*AO *CT A * ^ t RO tC 1 1 UM* VM ,CC » CH» 

frf- fCV.rp.O* ,PV,::*,CT,TAUV,TAUV,TTOT,N,0,MU ,T,TN,rr ,TAM,TAIR,0,C » 

M K.OC 

f.? CALL »’HPiPAICT,Ch,CV,CP,Cr,0X,0y,D2.0T,TAUX#TAUY ,T TO T ,OR ,ff ,RRf 

6! CHH,MV»?H,l?V,TAIO) 

t.U CALI TRY TAIT tJtlTit*iN>CrA ) 

t.*, CALI Pm'Vli.J.HtlKtJKfMNtHtC) 

t,t CALL r(^k'(IN,J*i,KN,WI 

1,7 CALL PRTE“I TN, JN,K\,TNI 

♦,« CALI PKINtPIIM,JN,KN,P| 

7.9 call P^I^TH(IAUJ^,NI) 

7'- STOP 

71 f NO 
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SUBROUTINE PROGRAMS FOR NASUM III 
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7. >.2. B0UND2 

This subroutln* sets th« boundary condition* for 
all varlablaa. The boundary conditions for velocity are no 
slip and no normal velocity on the bottom and on the shor<!;' 
line (y-axls) . For temperature, the bottom and shore line 
are treated as adiabatic. For oper boundaries, the boundary 
conditions on velocity and temperature are -0 and 

*0 respectively. The temperature boundary condition 
on the surface Is specified In temperature subroutine "TEM5". 
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6 
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B 

9 

in 

11 

12 
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IS 
1 

17 

18 
19 

22 

2Z 

2<* 

•> 

2 A 

27 

2^ 

2*^ 

in 

ii 

32 

33 

34 

35 

36 

37 
3*^ 

39 
4 C 

41 

42 

43 
4 4 

45 

46 
4 7 

4 A 

40 

50 

51 
5? 

5 3 

54 

55 
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SUBf^OUTINC B 0 UKn 2 ll«J«H*lN«JN»KN«U«V«H«C»D 9 r fWfHitMre tTfTN»TFI 
OIHfNSION U(1N vJN tKN) fVi INyJNfKN I fHIINfJN fKM t 9 
CGUNtJNtKM , 0 ( IMfJNtKNIf C 9 TIIN 9 JN 9 KNU 

CT^M TN»JI4«KNI fTTI I »Jh 9 HN I 9 Wll N, UN tHM I ,H I <IN 9 JN 1 » HTC I IN^JN I 
KNl =KN -1 

IM riN-l 
UNI ruN -1 
00 S ICI 9 HN 
00 S l:l 9 lN 
C 

C TAR U BOUNDARY 

C 

ni 9UN9K |:ni 9UNI 9KI 
TN( T9UN9K UTNdtUNltKl 

TFa 9 UN 9 K)=TF(l 9 UNUK) 

HTECI ,UN)=MTEd,UNU 
D(l 9 UN,Kl: 0 n 9 UNl 9 KI 
Lil 9 UN 9 K |:£ (1 9 UNI 9 KI 
W(I tUN^KlrUCItUNl ,KI 
5 CONTINUE 

C 

C FAR I BOUNDARY 

C 

00 10 Krl ,KN 
DO 10 U=l|UN 
r(IN,U,K)rTriNl,U 9 H) 

TM INrUfK irTNCINl 9 U 9 H I 
TrnN,UfK irTFUNl 9 U 9 KI 
HTE riN,U)rHTE(iM ,j| 

D ( I V , U 9 K ) =0 (I N 1 , U 9 K 1 
r. (INiU.K )=F i INUU 9 HI 
i^'dNfUvK )rwriNl 9 U 9 K) 

1C CONTINUE 

C 

C ALONG I BOUNDARY 

C 

DO 15 KritKM 
DC 15 I-1,IN 
Tn»l 9 K)rT(I, 2 fK) 

TN(r,l,K)rTN(r| 2 »KI 
TF( I 9 I tK iriF (T t2,K I 

H re ( 1 , 1 i=mtf(i, 2 i 
0 (I ,1 ,K ) ;: 0 ( I ,2 .K 1 

F(i,i,K>=!:n,?,K) 

W (I ,1 ,K irwn ,2,K 1 
15 CONTINUF 

c 

C ALONG SHORE LINE 

C 

DO ?0 K = 1 ,KN 
DO :n u=i,uN 
in ,U,K 1 =T( 2 ,J,K ) 

TN( 1 ,J,K )=TN( 2 ,J,K ) 

TFU ,U,K ) rrr (2,U,K 1 
MTE ( 1 ,U ) rMI ( 1 , Jl 
0 n , J,K ) = 0( 2 , J,K 1 


57 


C 11 *JtN 1 

58 


W(ltJ»Ki:w<2»J»KI 

59 

2P 

CONTIMUC 

60 

C 


61 

C 

BOTTOM BOUNDARY CONDITION 

6? 

C 




00 2S I=1«IN 

6a 


DO :S j:i,JN 

65 


T(I,J,S»=T(I»J**tl 

66 


TNI I»J»5i:TNlT«J»AI 

67 


TriI,J,S):TFII,J,RI 

68 

2S 

CONTINUE 

69 


RETURN 

70 


END 


■'Ht-CEOING PAGE BUNK NOT FILMED 
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7.2,3. CONST 

This is a small main program and need to be run 
In order to determine the constants DEEX, DEEY, EEEX, EEEY 
which are used In the subroutine "XYSH". The Input cards 
reads 


(1) 

XB, 

A. 

DX, 

AN 

(2) 

YB, 

B, 

DY, 

BN 


where XB, YB are X and Y boundary distances. A and B are 
distances from the shore line (y**axis) and x**axls respec** 
tlvely from where the stretching starts. It Is shown below. 

V 

1 


■ 







^ V 


DX, DY are the minimum grid size needed. AN, BN are the 
number of grid points in the x and y directions. 


X 

2 

3 

<1 

5 

6 
7 
A 
9 

10 

11 

12 

17 

IM 

15 

16 
17 
la 

19 

20 
21 
22 
23 
2<* 

25 

26 
27 

2a 

29 

3H 

31 

32 

33 

34 

35 

36 

37 
36 

39 

40 
4 I 

42 

43 

44 

45 

46 

47 
4ft 

49 

50 

51 

52 

53 

54 

55 

56 


50 RtAO 4,XBfA|0X9AN 
4 FORMA T( ) 

1F(XD«LF.0«I GO TO aqO 
yRITEf6»lll XBfAtDXvAN 

11 FORMAT ( • 1 •♦•XBK0RV=%F10, 1 ASSFlO.l *• 

INBR PniMTS=*#Fl0.1/l 
EUXR-A 
E2-IA^I-1.M0X 
OC=*Q 19A 

16 CrO* 
j:0 

8 CrctOC 

J — *1 ^ 1 

0::C«AL0G (A/C«SCRT HA/C l♦«2♦l • M 

UrAMiSH <E2-0>/Cf65,l 

X=A^C»5INH(UI 

WRITE«6,42I C,X 

IF(V,nT.XPI GO TO 8 

IF< J.r.T. n GO TO 17 

oc=nc/ 2 . 

GO TO 16 

17 CMIM=C-DC 
CMAX“C 

33 CUNTI^Ue 

Cll =C^AX 

0=C1 l>iiALOC(A/Cll •SOPTC U/Cll )♦♦2♦l.) I 
urA'MNl Hr.2-01 /Cl 1 ,85.1 
ERRU-CltCll»SINH(Ul 
J-1 

CK =CfAX-0C/2. 

1 Jrj4l 

IF (J. nr. 30) GO TO 99 
OrC12»>ALOO(A/C12^5 0RTHA/C12)**2^1.n 
UZA'-INI I (E2-C)/C12,85.) 

ERR3r-r 1 ♦C12«S INH(U) 

WRITE<6,42) C12,D,U,ERR2 
42 format f 1 X ,40 15.7) 

IF( 5 ( FRP2/XP ) .1 T..C01 ) CO TO 2 
C13--<cil«FRP2-C12*tPRl)/<rRR2-eRRl) 

C13rA^AV I (C)3,CMIN| 

CnrAMIMl (C13,CMAX) 

cure 12 

C12rCl 3 
EPR 1 rFRP2 

uo ro 1 

2 ClrC12 

hIRI TF <6 ,6 ) Cl ,0 

6 F0RMATHH0,*Cl.»,ri5.7,» Or*,Fl5.7/) 

nUM rn . 

NriNU Am 
00 3 lrl,N 
XLr(l-l )^0X 


DELTA X 



S7 X:A«Cl«SINHnXL-OI/Cll 

5P OCLTASX-OUM 

19 DUMSX 

ec' wRirei6»5i TtXLfXtOeiTA 

61 3 CONTINUC 

67 5 fORHATnXt*ts*»I*»#* XLS*,rjO«?#* XS»tn0.2»* DELTAS*, fio.l I 

63 GO TO SO 

6« UR1TEI6,1«) 

65 rORMAMlX.’NBR ITfRATICNS EXCfCPEO 30*1 

66 GO TO 400 

67 4.1C STOP 

6^= END 



'.’in'tn-DtfiQ PAGE BLANK NOT ^^L!V^LL ! 
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7.2.4. CURNT 

This subroutine sets the velocity field in the 
whole domain equal to the current velocity. In this case 
2 cm/sec is chosen. If the initial current is more, the 
values should be made equal to the measured value of current. 
If there is no initial current, this subroutine can be de~ 
leted from the main program. 


SUH^OUTINC CliRNT ( 1 »J »K , I N ,UN ,KN tU t V t H *6 tO *C ) 
OlHriSSlON Ul |N,«IN«KM »Vi IN,«IN«KNI 
r,H( IN, JN,KN) ,GIIN,JN,KM ,0(INtJN,KNI *C(IN,JN,KNI 
KM :NN-| 

(lO 10 K = 1,KN1 
UO 1C 1 = 1, IN 
00 1C J=1,JN 
U( I ,J,K I =C.O 
V(I ,J,H |:?.0 
Nil ,J,K I =fl.O 
MI ,J,K 1=2.0 
Dll «JfK 1*C«0 

CONTlNUf 

Kf.TURN 

LND 


original 

POOfi 


PAGE 


IS 


^oality 
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7.2.5. DENSTY 

This subroutine computes for the whole domain the 
density using the temperatures computed In the subroutine 
"TEM5". 
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I 
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1 
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1C 


SUDnoUTINf nCNSTVI lN,JNtKN,POtT) 

OIMTSSICN ftO{TNt>IN*KN} *T (TNfJNyKNl 
CO 10 l:l,IN 
00 1C JsltJN 

00 to KritKN 

wot t *JtX ):l.0n092P-0*0C0C19*IT(I»JtM I )>0* 0000096 * 
rmT».j*t()*Tii.j»Kt) 

1C CONflKUE 
WCTURN 
CWO 


1 

1 

1 


ORIGINAL PAGE IS 
OF POOR QUALITY 
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7.2.6. ETT 

This subroutine calculates the wave height by 
subtracting mean level of water (HI) from the total height 
computed. 


n ^ 


THIS PROCPAH CO*4AUT(;S THE UAVC'MCIGHT 
SUSROUTIMC CTTITM*JN*HT,HI,HAr,FTn 

01 HENS! ON HTI1N*JNI fHl ( IN t JN ) ,HAR ( I N > *C TA ( IN* JN) 

00 10 I:i*IN 
00 10 JSliJN 
CTAfItJ l=HT(I *JI-H1(I,JI 
10 CONTINUE 
HE TURN 
ENO 


iEOEDIKG PAGE BtAW! NOT FILMED 
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7.2.7, H1TE2 

This subroucine computss th« fret surface helghc 
using the equation 


/ 

O 


'3h 

*b * “b ^ 




The nvttaerical scheme used is forward difference in time 
and central difference in space (FTCS) . Simpsons rule is 
used for numerical integration. 
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IP 
1^ 
:o 

2 \ 

n 
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Z I 

:<i 

2^ 

2b 

27 

2«» 

2« 
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SI)srV0l)7 INC HI TC2 < IN,JN »K N ,NAI? tU t V,NT «HTO»HTr tHZtOT tOX »0V ,HOUN 
C»XX , VV tWN I 

ION N«HI IN, JNi »U( INiUNiKN I »V(TH ,JN ,NN I ,HT I IN ,JNI,HrOIIN ,JN I t 
CHrU>4| IN,^M •HTCI IN,JNI ,X X ( IN 1 .V V ( JN I ,Whl IN, JN,NN I 
KNH 1 =MN-1 
IMl rIN-I 
JM rjN- 1 
DO 5C I:2,IN1 
DO SO J72,JNl 
HDUH(I,j|:0.0 

00 60 h;i,kn 

IF<^A<»IIiJ>,eC.Ill GO TO II 
..0 TO 50 

II niHUX: (MT (1 «l , JI«UII* I ,J ,K |*HT( I2.«0X I 
01HVV:|HT ( I • J*n *V<T, J* 1 ,HI*HT( I , J*l I *V( I I *M I 1/ (2 . «0r I 

C.... SINPSCN'S suit IS usro FOR INTEGRATION 
IF(K.C0.1 .CP.K.CC.S) GO TO lOI 

iFiK.'-o.r.cp.x.ro.si co to ios 

r ,j|;( (0 |HUX»XX( I MO 1HVV*VY<JM«0Z*I2./3.M •H(XIMlI,4t 
00 TO 103 

101 MPU“( I,J»rnOtHUX*XX(I M0IHVV*YY(jn*02/3.MNCUM(I,«n 

00 ro 103 

1C2 I ,JM( (01HUX*XX( I M01HVY«YY(JIM02*M./3.)I «HOUN(I ,JI 

If 3 CONTIVUC 

HU n , j» = MTn « r ,j) -HOUMd , jmot 

6C: CONTIMUC 

HTt (I .JMHTE ( T,JMWH( 1 ,J,KNMOT 
SO coNTi*iur 

^rrt'RN 
tNO 


7.2.8. INITIB 

i*hia Bubroutlne r«ads th« depths for a conatant 
depth baa In and Initlallzea the valuea u,v,w,p,p,7 and HI. 
The program sets u,v,w and wave height ETA equal to zero. 
The temperature la aet equal to the ambient temperature. 
The preaaure la hjrdroatatlc. 
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10 
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1«3 
2 C 
21 
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24 

25 

26 
27 
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3n 
!1 
3? 
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SUBPOUT INC TNITIDIINtJNyKNtC TA,H7«HI lUvVfPOfPfCP iRPtO^tHTOfO » 
CC.HtCfHTC *TtTN,TF tTAHf UH ,VM»W,Oir fOY I 
DIMCNSION £TAIIN,JN) ,HT( I*/«JNItHll|N9 JNI «U( TNtJNfKNI vVdNtJNvKN) • 
CPOl TNtJNtKM t^flNvJNtKNI ,H7n HN vjM tW (INtJN vKN) 

CiDllNt JN tHN) yfllM »JN,KfJ. , H (I N 9JN 9KN I • GUN tJN 9 K N) fHTf ( IN» JNI 
C»TAMnN«JN,KN) fT IXN«jN|KN)»TN(Xr4,JN»KNUTFnN,JN»KN) 

C«UM (IMtJNtKNl tVHUNtJNfKNI 
00 500 JSliJN 
HI U »J1 :3C0«(1 
HH ?t Jl=350«0 
HII 3fJI=40a»0 
HI(4 tJ)::4S090 
Hl(5,JlrS0C«0 
HI(6,J)r550.0 
Hn7,J)=6CC.O 
HI(«,J|r6C0.0 
Hl(9,J|r6CO.O 
HM ICf J)::6an«n 
Hl( n ,JU600.0 
HI( !2 f J)=650,0 
HT( 13fJ)r700.0 
HI( l4tJ)=750.0 
HM l5,JlreOO.O 
HI( lb,JJ:R50.0 . • 

HU J7,J)=900.0 
Hli \8 9JJ r95n.n 
Hi( ipfjiziooa.o 

Hi(:otJ)=iooo.o 

500 CCNTINUF 
00 in 

00 1C J=l»JN 
ITA { T ,jjrc.n 

HTD< I I,J) 

HTC f 1 ,J) =HI ( I ,J) 

1 r CON r iN'UF 

DO P iri ,1^ 

00 »» Jll , JN 
DO K = l,KH 
UM( T,J,K)r0.n 
VH( T , J,K ) =0.0 
w n ,j ,K ) =G.n 
U ri ,J ,K I rUM ( I , J,K 1 
V a , J .K J =vf;c I . j,K I 

TAM ( I , J,K ) = 2 5.C 
T (I , J.K 1 =TAM ( I ) 

P ( I »J ,K ) =CP*HT (I , J l♦PR♦< K-n 902 
H (I «JtK ) =U( I « JfK ) 

0 n » J »K ) =tM r t J »K ) 
fi ( I , Jf K ) =V( I , J «K 1 
£ (I ,J.K ) =V( r , J »K ) 

T N ( r , J , K ) I T ♦ I , J , K 1 
TF a ,J|K ) rT ( I I 

« continue 

RETURN 

END 


•rli'ir^AL iPAGfe i\6 

POOR QUALITY 
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7.2.9. INLET 

This subroutine puts the velocity and temperature 
of the discharge at the discharge location. The value of 
velocity specif -'ed in this subroutine should be calculated 
depending upon the mass of the discharge. 


SUnnOUriNE INLET(I,J,K »IM,JN»KN,WH,T,TN|Tn 
OIMEN', lOM r(lM,JN«HM tWHIIN,JN,KNI 
C«TN(P.',JN.KNI ,TF( IN,JN,NN) 

DO 10 I=8»10 
00 10 J=in,12 
00 10 KrlfKN 
VHII.JtSir-O.JS 
T (1 ,J,K 1=35.0 
TN( I ,J,H 1=35.0 
TF< T,J,K 1=35.0 

continue 

HETURN 

END 
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7.2.10. MIXT 

This subroutine mixes the temperatures by an aver- 
aging process in such that unstable density gradients are 
eliminated. 



I 

I i 

4 

I ‘r 

0 

9 

10 

I 11 

1? 

14 

15 
14 
17 
n 

19 

2C 

21 

2 

2 " 

24 

2^ 

2t 

21 

»? 
3'^ 
31 
3 ^ 
3 3 
3 4 


sun POUT I Mc Hixrn 9 J»h »iN»jNfHN«n 
0IMC^<5I0N UlNtJNjKN) 

00 !C lrl,lN 

00 10 jrl^JN 
Kri 

If <Tn,j,KKGc.Tnfj*K«in go to i 

AVTr « T n tJfH MT ( r «JtK« 1 n ;2«0 
HI 9J,K)r^VT 
T (I tU»K#l > = AVT 

1 CONTINUr 

ir nn»J,K#l).GE.T(l»JfK«21 I co to 2 
AVr = (T(ItJfK)«TfI,J«K«l)«nr,J«K42l)/3«C 
n I fjf K ) rAVT 
T (I , J,K^ I UA VT 
ni )rA VT 

2 CONTINUE 

ir ( T n ,J tK ♦2) .GE .T( I # J,K*3) ) GO TO 3 
AVfr<T(T,J,KMTn,J,K^lMT<ItJf‘<«2)^Ta,J,K^3n/4,0 
MI ,J,IO r^vT 
M I I J » K ♦ I ) r A V T 
MI ,J,K ♦: MA VT 

1 (I ,J,K43MAVT 

i COMINUr 

Ir <M r »JfKt3l .C£ .T(l tJ#K^4n GO TO 4 

AVT = ( M r , J,K ) ♦ T< I , n ♦MI ,JtK ♦? I 4T (I »J ,K ♦! ) #T (I ,J,K#4 n/S.O 
T n , J .K MAvr 
MI ,J,K ♦ i )-AVT 
M M J 9 K M' ) = A V T 
TM ,JtK ♦3UAVT 
M I 9 J, K ^4 UA VT 
4 CuNTlNUE 

r' CONTTNU^^* 

RETURN 

END 
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7.2.11. OLDHT 

This subroutine sets the values of height HTD at 
time level n to HT at time level n-1 and HTE at time level 
n+1 to HID at time level n after all computation are comple. 
ted. This Is necessary In order to retain values of height 
at one time step lag. 
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1 

2 
J 
U 
5 

(i 

T 

a 

9 

to 

ii 


C THIS rROCRAM TPAHSFCRS MATRIX MTO 10 MT 

C 

SUMPOUTINf OLOHTnNfJN«HTE fHTDtHT I 
OlhfNS ION HTD (I rj,^N ) rMfE ( IM »JN) tHT (IN»JN) 
00 lU 

00 10 J=1«JV 
HTC If Jl =HTD( I f Jl 
HTO(I» JI^HTCdtJ) 

10 CONTINUE 

RETURN 
END 
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7.2.12. OLDUVT 

This subroutine sets the values of velocity 
and temperature at time level n*fl equal to the values at 
time level n and the values at time level n are made equal 
to the values at time level n*l. This is necessary in order 
to retain the values of velocity and temperature at the 
time step lag. 
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t 

2 

5 

H 

5 

6 
r 
a 

9 
13 
11 
12 
1 3 
l<» 

15 

16 
ir 

IB 


i 


C THIS PR06BAH TRANSreOS MArKTCCS U*V TO 0*C RfCSPCCTIVCLV 

C AND T TO TN 
C 

SUBROUTINC OLOUVT (1N,JN tU(V,H,CtOtf tT«TN,Tri 
OIHCNSION Ul 1 Nf JM«KN l,V( 1N,.JN,KNI ,0«lNyJNtKN) rC I IN,JN,KN) , 
CH<INt>iNtKNI«R(lMt Jri,(;NltT(lW(JNtKNI|TH(Ii4iiJN»KN) , Tr ( 1 JN ,KN I 
00 to K:1«KN 
00 10 l=lflN 
00 10 UrltJN 

V( ItJfK) -C (I f J t K ) 

H( I. J,K)=0< 1 1 J*K I 
GI1«J»H)=EII *JfK ) 

TU«J»K i:TN< I *J,K I 
TN( I, J«K)=rFn, J,K| 

10 CONTINUE 
RETURN 
END 
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7.2.13. PRPARA 

This subroutine prints the input values and the 
total time the model is simulated. 
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I 


1 

2 
3 
H 

5 

6 

7 

8 
<1 

ICi 

11 

12 
n 


C THIS SUOROUTIM: prints P8RAMCTCRS FOR THC FRCF SURFACt MODEL 

c 

SUanOUTINC PNPAf^A (CI »CM|CV»CP »CC tOX,OYtOZ*or,r«UXtTAUy| TTOTfGRtPr t 
CRP,KHtKV,PH*nVt r«IR) 

PPINT I»CI,:h,CV,CC,CP,OX,OY,DZ,DT,TAUI(,TAUV,TTOT,G»,FF,RR,KH,<V 

C,BM,6V,TA1» 

I FORHATt/* Cls»,n5.7,/' CH = • ,C 1 B . 7 , /• CVr»,n5.7/» CP = »,EJ5«7, 

C/» CC=*,E15.7,/» DX = *,eiS,7,/» OY = • ,C IS .7 , /• DZ= • ,t 15 .7 , / • OTr% 
CE1B.7*/* TAUX=»,tl5.7,/* TAUY=* ,CI5.7 ,/• TTOTr • , [ 15.7 »/ • Gk=», 
CE15.7,/* FFs* ,E 15.7t/* RR = • »E 1 5 . 7 ,/ • KMS», E15.7,/» KVS»,EIS.7,/ 

C* BHZSE1S.7,/* BV=» ,ri5.7,/* TAlR:*,Fi 5.7/1 
RETURN 
END 




•y 

/ 
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7.2.14. PRES 

This subroutine calculates the pressure field 
using the updated density field. 
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1 

2 

3 

1 

5 

6 
7 
1 
0 

1C 

11 

12 
13 


C THIS PROCRAH CAICULATCS THE PRESSURE PICLO 

C 

SIIPROUTIHE PRPSIlH»JK*KH,HTtROtCRtP«0?l 
OIMENSION HT IIN,JHI,tiO(lN»JN»KNI ,P(lN,JNyNNI 
00 10 IriylH 
00 10 jriyjN 

pix yj*! iso.n 

00 8 MS?yAN 

P(l yJyN ISPIT yJyN-l l«6R*HTt IyJI«IROI I yJyN-ll«RO(I yJyK I l«PZ/2»0 
a CONTINUE 

i(. continue 

RETURN 

END 
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7.2.15. PRETA 

This subroucine prints wave height (ETA) . 


u o 


THIS PRObRAH PRINTS THC WAVE HEIGHT 

1 SUbROUTINC PRETA(ItJ(IN»JN»CTA) 

<1 OIHENSION ETAUN.JN) 

5 DO 10 I=ltIN 

S 10 PRINT Il•lt<ETA(l«J>,J=lrJ^<l 

7 11 FORMAT*/* 1 =’,I3/» MAVE-HEIGHT*/(5X,8EIS.711 

B RETURN 

9 END 
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7.2.16. PRTEM 

This subroutine prints temperatures in the whole 

domain. 


3" VI 
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7.2.17. PRUV 

This subroutine prints u and v velocity in the 


whole domain. 


1 

2 
1 
4 

s 

b 

7 

f) 

9 

10 
1 I 
12 
1 ^ 


C THIS PROGRAM PRINTS TMC HORIZONTAL VELOCITIITS 

C 

SUUROUTINC PPUVII tJfKtlN tJN»KN«UfV) 

UTPrNSION UdNiJNiKNl fVnN ,JNfKN) 

Kra=KN-i 
UO 10 KSltKN 
UO 10 iritiN 

PRINT cud «JNf 

10 PRINT 12f(V(ltJ«K UJzitJNI 

11 FORMAT!/* Kr*,X3,3Xt*I=’*I3t/* U-VCLOCITY • / C 5X t8E15. 71) 

12 format ! * V-VELOCITV*/ (SX t8C15.7) ) 

RETURN 

LNO 
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7.2.18. PRW 

This subroutine prints n values in the whole domain. 



7.2.19. 
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pm 

This subroutine prints w-veloclty In the whole 


domain. 


suePOuriNC prwh(in ,jn,kn,whi 

OTHfNSI^N WH(IN,JN,KN) 

DO 10 Kri,KN 
00 10 I=ltlN 

PRINT 1 ItKtIt (UH( I» J»K) ,jri,JNI 

FORPAT</« K=*,I3|3X#» I r*,I3t/* NH-VELOCITY*/ (5 X , BFl 5 . 7 ) I 
RETURN 
END 
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7.2.20. READl 

Thlfl subroutine Is used to read the Information 
stored by subroutine ’’SIORE". This subroutine Is used 
from second run onwards In order to use the values created 
in the first run (ie IRUN-0) . "READl" and "STORE" corres- 
pond to each other. The "READ I" subroutine uses a file 
designated as "UNIT 7" In order to read the Information 
on the tape. 
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I 

3 

4 

6 

7 

0 

9 

10 

11 

12 

13 

14 

15 

16 
17 
IH 

19 

20 
21 
12 

23 

24 

25 

26 
27 
2H 

29 

30 

31 

32 


C this SUHROuTINC RE:a05 DMA TROM TAPE 

C 

bU( ROUTINE RCAOl (IN «UfV«W«UI tHT fHTOfMAR vFTA|R«RO»CI*U*^» VM» 

CCC»CH,C Vf CPiDXf r)Yf02frT, TAUX,TAUY vTTOT tM,G»HTE ft ttUfTF, TftM, TAI^tOt 
CCI 

DIKE NS ION iJ(lK»JN,KH ItVflNtJN^KNI |M (I N » JN »KN 1 »P ( I N, JNf »^N » t 
CMK IN, UNI tMT C IN, JM tHTD UN tJN I ,NAR i IN tJN I , 

CETAI in, JNI,R0< INf JNfKNI ,HUN, JN,KNI tUllN, JN^KNI ,HTEI in, JNIt 
CTUN,UN,KNI, TNaN,JN,KN I ,Tr (IN, JN,kN I ,T ANC IN, UN, KNI 
Cf UH( IN, JN,KNI ,VN (IN,JN,KN l,n(IN,JN,KNI,EnN,JN,KNI 
1 CONTINUE 

READ! 7,EK0 = n ( ( UM I , J ,K I ,K r 1 , KN I , j: 1 , JN I , Ir 1 , 1 Nl , 

Ci ( I V( I ,U,K|,K :l ,KN| ,Url ,UN I ,171 tXNI , 

C( M W( T,J,KI ,K=1,KN| ,jrl ,UN 1,1 71 ,INI , 

Cl UH( r, J,KI ,K 71 ,kNI ,jrl , Jn I ,I71 ,INI , 

C( ((G( I , J,KI ,H 71,KNI ,J71 , JN l,lrl ,INI , 

C( UD( T, J,H),K7l ,KM ,J7l ,UN I,I71 ,INI , 

C( i (£ ( fINIy 

C( ( (‘MT. J.K»,Hrl,Kro »Jrl,.tM,I=).lN). 

C( ( <r;u( 1, J,K> rK = l »KN ) ,J = I , JN I tl = l t 
C( ( ((iM( 1 , JtK> tK = l tHN ) ,jri , JN = 1 , IN), 

C( ( « VM( 1 1 J,KI ,K:I ,KN ) , jri , JN) , izi , IN), 

C( INTO < 1 ,J) ,Jzi,JN ),1 ri ,IN ) , ( (HIE ( I ,J) ,jri , JN) , ir 1 , IN) , 

C( (Hid, J),jri,jM) tlrl.IN) , l(HA9(I,J), J=l. JN),lri ,IN), 

C( (HT( I, J),jrl ,JN) ,I ri ,IN ) , ( (ET A ( I ,J), jri , JN) ,lzl ,IM , 
C((mi,J,K),K=l,KN),Jzl,JN),I=l,IN), 

C( ( ( TN( I, J,K) ,KZ1 ,HN ) ,J=1 ,JN ),IZ1 ,IN), 

C( ( I TF( I, J,K) ,K = |,KV ) ,JZI,JN),I=I , IN), 

C( ( (TANd, J,K ) ,K-1 ,KN I ,J = 1 ,JN) ,iri ,IN) , 

CC I ,CC,CH,CV,CP,OX ,Or ,Ojr,OT,TAUX ,T AUY,TTOT ,TAIP 

(RETURN 

END 
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k 


7.2.21. READ2 

This subroutine reads the "MAR" matrix. The MAR 
numbering system’ Is used In order to Identify the points 
In the interior, on the boundaries and outside the domain. 

The "MAR" numbering used is as follows, 

MAR (I,J) * 0 for points outside the domain. 

MAR (I,J) - 1 for upper horizontal boundary. 

MAR (I,J) "• 2 for lower horizontal boundary. 

MAR (I,J) - 3 for left vertical boundary. 

MAR (I,J) - 4 for right vertical boundary. 

MAR (I,J) - 5 through MAR (I,J) - 10 are boundary corners 
and are specified as below. 
r-lAR (I,J) ■ 11 for points in the iterior. 

MAR- 5 MAR-6 MAR-7 MAR-8 MAK-& MAR- 1C 

P L L. ”1 I “I 
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I C IHIS &Of<40uriNC Dcrires the HAR MARTRIX for L0CATIN6 THE POSITIONS 

? C 

1 SUl<'70UTTNE PE«0;nN*JN,HAR| 

R Q1H<'N1^I0N HAIM INfUNI 

5 MARIlit):? 

6 HAn«l,JN|:S 

7 MAiMIWfUSR 

8 HAMIlN,JN)StO 

9 1 NM I S I N - J 
10 

11 00 13 I=?«IRr<l 

12 *'APir,n=2 

13 1(1 «‘AR«r,JN»Sl 

H PO 2C Jr2,jNM| 

15 maRUiJ)::3 

16 2C PAPn*j,JI=H 

17 00 (0 Irr.INHl 

18 00 ?u 

19 I" MAR n , J > =1 1 

20 PrTI'RK 

21 LNP 
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7.2.22. STORE 

This 8ubrouCln« Is ussd to score the valuitB 
at the end of all computations on a file designated as 
."UNIT 8". 


INIS SUBROUTINC STURCS THE DATA INTO TAPE 


SUIlROUTINC STORE (IN ,JN «KN »Uf V fK »HI .NT tNTO »H AR ,£T A «P , RO, C 1 1 U H, VH, 
CCCtCH»CVtCP»OX»DV«DZtDT •TAUX,TAUV»TTOT»H»G*HTE»T ,IN,TE» TAM, TAlRf 
CO*EI 

OIHENSTON UUN,JNtKN I.V(IN»JN,KNI »W (1 N, JN .KN I ,P ( IN»JN»KN) , 

CHT( IN, jrO.Hl (1N,JS) ,HT0(1N ,JN l,^<AP (IN.JNI ,0(IN,JN,KN| ( IN, JN,<NI • 
CETAdN, JNI,R0IIN,JN,KN»,HI1N,JN,KN),6I1N, JN,KN|,HTE(1N, JN), 

CT( IN, JN,KN>, TNdNtJN ,K N ) , TE (IN , JN ,K N) ,T AN (IN , JN , KN ) 

C,UN( IN, JN,KN» ,VH ( IN , JN,KN ) 
hPITE (8H((U(I,J,KI,K = 1,KNI,J = I ,JN),I=1,IN), 
C(((V(I,J,K),K=1,KN),J=1,JN|,I=1,IN), 
C(((U(I,J,K),K=l,KN),J=l,JN),lr|,IN>, 
C(((H(T,J,K),K=1,KN),J=1,JN),I=1,IN), 

C( ( (G( I ,J,K) ,K*1, H^i) ,J*1 ,JN ) ,1*1 , IN I , 

C( ( ((HI • J,K) ,K = l,Kri) ,J = l, JN I ,1=1 ,INI , 

C(((E(I,J,K),K=1,KN) ,J = l,JN > ,1=1 ,IN> , 
C(((P(I,J,K),K=1,KN),J=1,JN),I=1,IN|, 

C( ( (PO( I, J,K) ,K = 1 ,KM ,J = 1 , JN ) ,1 = 1 , INI, 
C(((UN(I,J,K),K=1,KN),J=1,JN),I=1,INI, 

C(((VM( I, J,K) ,K=I,KN) ,J=1, JN),I=1 ,IN), 

C( (HT0( I, Jl, J = 1,JM, 1=1 , IN ) • ( (HTE( I,J) , J = 1 , JN I , 1= I , IN ) , 

C( (HI( I, ,J=1 ,JN I ,I =1 ,IN > , ( (MAR ( I ,J),J=1 , JN) ,1=1 ,JN) , 

C( (HT( I , J) ,J=1 ,JN) ,I=l,IN),((ETA(I,J|,jrl,JNI,I=l,IN), 

C( ( (T( I, J,K) ,K=I ,KN > , J=1 ,JN| ,1=1 ,IN| , 

C( ( (TN( I ,J,K) ,K = 1,KN I ,J = 1 ,JN 1 ,1 = 1 ,IM, 

C( ( ( TF( I, J,K) ,K = 1 ,KN ) ,J = 1 , JN) ,1 = ] , INI, 

C( ( (TAM( I ,J,K I ,K=1 ,(<N ) ,J=1 ,JN) ,1=1 ,INI , 

CCI,CC,CM,CV,CP,OX,DT,n2,nT ,TAUX ,T AUV,TT0T ,TAIR 
END FILE a 
END FILE 8 
RETURN 
END 
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7.2.23. TEM5 

This subroutine computes temperature only In the 
Interior of the domain. The schemes used are forward In 
time and central In space (F.T.C.S.). Vertical diffusion 
term Is treated by DuFort~Frankel scheme. 
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i 

u 

c 

6 

; 

Q 

ir 
II 
1.^ 
I ? 
1 

IS 

u 

1 7 

2: 

22 

2*^ 


7 ^ 
27 

2^ 

J 

31 

3.7 

33 

5^ 

T*' 

37 

3 ^ 
3^ 
u r 

M i 

4 7 
4 ^ 

4 4 
4S 
4f 
47 
4^ 
(40 

r.o 
s i 

5 2 

53 

54 
s? 

56 

57 


fU^ROUT IM: TrM5 < IN,JN,. N,HT,Hfn.MTr fDX , 0 Vtn 2 , 0 T tPHiBViTtlNf rr, 
Cb »H ,G,MAR ,HK ,TAIR ,TAM tXX,YVtXXXfYyVtLtlNI 
DIM^NMON HTUNtJ^ltMTOf lN«JN>tHTrnN«JNi •HlPdNtJNIt 
CT^MflNt JN»K^f I • b (XNvJNfXM } ,H« IN, JN^KN I tC IJN, JN 9 KN I ^RO ( IN « JN«KN 1 1 
(XX( !M , YV ) ,XyX (IN) , YYY ( JN I , 

CY n^', JS,KN) ,TN(INf JN,KN) ,TF( INtJNfXN ) 

K\l TKN -1 
1 M 7 1 ^ - I 
JM =J^-l 
00 ir. K 31 ,KN 1 

00 1C I3?,1N1 
i)C in J32,JN1 

rF(v/\p( i,j),c:o.ii) CO TO 11 

GO TO 3?C 
1 ! CONTI vur 

OHUTXr ( MTO( I ♦ 1 » J I ♦H( I ♦ 1 , J ,K ) tTN ( !♦ I , J ,K l-MTO CI-1 ,J) 

C^HU - I , J,M )<iTN (I -1 ,J,K n / (2^0X> 

OMvTVr CHTOn , jtl l^Cdf I ,KI bTN( I,J# 1 fK !• 

CH TCM I ,J- n^r ( I ,J-1 ,M ) »TN II ,J-1 ,K I ) / (2#0Y) 

0 ) vriHT II ♦J^1)-HT(I,J-1) )/l.?«OYI 

OIT Vr I TNI I, J4 1 ,K ) -TNI r, J»1 ,K ) )/ (2*DY 1 

0:TV=lTNII,j4l *K ) ♦TNir ,J-1,K )-2tTN(I ,J,K) )/f DY»OY) 

OHX = IMT ( I ♦! ,J) -HT I I-l , Jl) /(2*rx ) 

OIT vr I T^ I !♦ I ,J ,Kl -TNI I -I ,J,K n * (2*0X 1 

0?TXrc TN ( !♦ 1 ,J ,K ) ♦ TNI I -1 ,J,X )-2^TN(I ,J,K)I / (OXtOX I 

lFlM.ro. 1 ) no TO SC 

OIWT JTIW II ,J,K ♦IJbTNII ,J,K4l) -Wl )*TN( I,J,K-J I )/C?tOZ) 
DlT*’r|TN(T,J,K4l )-TNCI,J,K-l I )/l2*07) 

02T d 7 I TNI I , J,K4 I ) 4TM I , J,K-1 l-TII ,J*K ))/ I07«02) 

GO TO 2^0 

ir 01WT2 7 1 4«W| T,J,K4|)<»T^MT,J,K♦l)-5♦W(I,J,K)♦TNII,J,K)- 
ry 1 1 ,JtK 42 )«TM I , J ,K#2 ) ) / I 2^0 2) 

01 T23HT I I I TNII, J, n-TATR) 

orr .-’ii I ?«TM I , j,K 41 J-T (I ,j,K ) ) / lorv'-nz )- 2 *oiTZ/oz 

J' ; COST iNur 

LiMT7|MTr I T*J)-HTO(I*J) )/lOT) 

TLXrOHUTXbXX II ) 

TLY rn^.vrv^YY IJ) 

TL/-f<rD 1 1 tj ) *r \yj7 
lUdfK-l)«n7-n^PMTbniTZ 
Tl -:TLX4TiyML74TLC 

P' A V H I ,,HXbXX I I ) ♦CITY bXX ( I) 4MT I I ,J) ♦XX in*>XX C I 1 ♦D2TX4HT I I, J)<» 

rxx X 1 1 ) *n 1 TX ) 

II»Y f rtMY<rYY IJ ) ♦01 T Y»YY (J) 4HT ( I tJ) ♦YY IJ l^^YY IJ) ^nZTY^HT I I ,J)* 
TY YY IJ ) vP I TY ) 

Ij 0 TO s p 

S ')u I A* Tr X ♦ T!‘ Y 

TPVr|>VbP2T:i /MTO ? I ,J) 

T»'r Tfi ♦ TfJ V 

ir I T ,j,K I = 1 1 TR-TL T^nr^HTr (I , j)«tm I, j,K n / 

Cl MTr I T • J ) ^(<v^r T / 1 n 2*07 in 1 1 , J) i ) 

Of) TO fl 

300 TF I I ,J,K >rT^ M I I , J,K) 
o CCN r INUF 

crsTiNur 

RCrtlR U 
t M* 
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7.2.24. TIDE 

This sets the values of velocity at one boundary 
equal to the value of current comming into the domain. 


SUt! ROUT INC TIDE ( 1 *J,K , 1 N t JN,KN, U , V»H ,G *0 tt *T »TN« TF I 
dimension U(IN,JNtKNI «V( INyJNtKNI ,K( IN,JN,HNI 
DIMENSION C( IN«JN»MN| • 0 t I N » JN »K N I (I M , JN »K N > 
DIMENSION T I IN,JN ,KN) , 7 N (IN,JN»KNI ,TF (lN,JNtKN) 

MNl TKN-I 

00 1 C I=ltIN 

DO 10 K = UKNI 

V(I,UKt:;.0 

6(1 ,1 ,K )r2.0 

C(I ,1 ,K 1=2.0 

CONTINUE 

RETURN 

END 


original 

OR POOR 


page is 

QUAUjy 
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7.2.25. UVELS 

This subroutine computes u and v velocities from 
the two horizontal momentum equations. The schemes used 
are forward in time and central In space. DuFort- Frank el 
scheme Is used on the vertical viscous terms. This subroutine 
computes velocities only In the Interior. 
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00 10 J=2,JNl 

no ^ Krl,KNl 

IFI-a; » I ,J» .EO.l 1 » CO TO 11 
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DHVriMf ( 1 ,J* 5 » -HT 1 1 ,J-1 ) I /<2*0Y I 
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7.2.26. WHTOW 

This subroutine calculates the value of W(^) which 
Is used In the model from the value of WH (w) specified. 


V 


SUBROUTTNf UHTOU(IN.JN«KN»MTO,HTC,HT « CT A , Ot IT ,W *W H« OX tOY *07*0 T tNAR« 
CXX, YY I 

OIMCNSION MT (IN,JNI,HTO( IN,JN)»HTE UN ,JN) »HAR(IN (JNI *01 IN,JN,KNI » 
CC <INtJN*KNI ,U(IN , JN«KNI,bH(IN»JN,KNI , t T A ( IN , JN ) « XX I IN i* Y Y IJN | 

KNl =KN-I 
iNlrlN-I 
JNI =JN-t 

00 10 k;2,kn 
00 10 isZtlNl 
00 10 J72»JN1 

1F(MAR( 1 ,JI .eo.l 1 1 CO TO «« 

GO TO 10 

OMX=«HTr«I*l ,J»-HTEn-l*JI»/l2*0X» 

OHVriMTf «I,JM)-NTtll,J-lU/«2*0Yl 
CTAXr UTAH* 1 ,J|-ETA«I-1 ,J)I/ 12*0X1 
C TAYriCTA (1 • J* 11 -eiA( I tJ'l n /f2*0Y) 

OHTr|MT'>« I, Jl-HTE ( I,Ji ) / lOTI 

W <1 ,J ,K t : (WN (I ,J*K )•( ( (K -1 )«02>1 l*0HT-0(I *J*K l*fTAX*XX( II 
C-E ( ’ f J|K )*CT Ay*VV I Jl«< h-1 1*02*01 It Jt K >*nHX*XXC I) *<K-l )*02* 

CE < I , J,K I *rHV*VY( Jl 1 »JI 
CON TTNUf 
KETURN 
E NO 
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7.2.27. WVEL 


This subroutine computes 
(^) by using the equation 



the vertical velocity 

+ ( x'2^ 

O 

V* ^ h T» V* ^h \ 


the numerical scheme used Is central In space and trapezoidal 
rule is for numerical integration. 
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C THIS PROCf^AH CALCUtATCS THC VCRTICAL VELOCITY IN THE GANNA COOROS 

C 

SUBPOUriNC WE LI IHvJN^KN f UtV t V| NT fOY 9DY t 07 f H AR t VX « YV tUH I 
DIMENSION ill IN.JNthNI tVllNtJNvKNIvH? f IN t«INI 9 U I IN 9 JN^KNI » NAR IIN9 JNI 
C9XX IIMI ,YYIJNI 9VHI IN9JN9KNI 
KNITHN -1 
INiriN-l 
JNl rjN -1 
00 10 lr 29 lNI 
00 10 J: 29 UN| 

OUHtO. 

00 9 K=l9KN 

IF(M|P( 7 , J| .£0.1 1 1 GO to 20 
GO TO IP 

2( OlHUXr |MTII« I 9J)*Ufl4 1 tJfK I -MTI 1-19^1*011-1 9 Ji HI 1/ f2«0X) 

OlHVYriMT II 9 J«1 l«V II9J«1 fK l-HTI I 9 U-I 1*VII 9 J-I 9 KI 1/ I2*0Y1 
ir (K.EO.ll CO TO 17 
iriHtfO.KM GO TO 17 

0UM=DUM#O2’i»IOlHUX*AXII MO 1MVY*YY (Jll /MTCI 9 JI 
GO TO 9 

17 DUM = DJ»>i ♦07*10 I MUX •XXII M01HVY*Y YIJI) / I 2*HTIT 9 JI 1 
9 CONTINUE 

OUmtOUM#W (I 9 UfKN| 

WUOr'J. 

00 « KrjjKNl 

IFMA^I T 9 J> .EO.l I I GO TO 40 

uO ro ic 

4 0 OIHUX IrlHTIMMJMUIMMJfK-l MMTCI-1 9 JMUII-I 1 1/ I2*0X I 
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C IMVYlr IHT I I 9 J^l MVII 9 J* I fX-l MMTII 9 J-l MVn ,J-l 9 X -1 I 1/I2*0Y I 
WUOrwUO^O?* 10 IHUX*XXI I M 0 1 IIU X 1* X X I I) ♦ 01 WV Y* Y Y IJl ^0 IH V Y1 * Y Y CJ 1) / 
<MMTii,jn 

-U 9 UfKls-i»U 0 ^DUH*(K-l 1*0 7 

a continue 

CONTINUr 

return 

CNO 
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7.2.28. XYSH 

This subroutine does the stretching in both hori- 
zontal directions and determines the constants XX(I), YY(J), 
XXX(I) and YYY(J). This subroutine needs the values of 
DEEX, DEE7, EEEX which are read in>the main program. In 
order to obtain these values, another main program "CONST" 
has to be run. 
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uo ic; UI,1*J 

XU l?U-l )4PCLX 

Ml )ztxi I i-f)irx)/rirx 
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DO ?C jri ,JM 
YIJISIJ-I l•^CLY 
BIJiriYfJI-nc^VWECE/ 

Yy(JUU/C05H<9( jn 
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Pig. 5 Grid configuration and MAR node values 
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Fig. 8 Horizontal Grid Point System Without Stretching 
For Free Surface Near Field Model Applied to 
Hutchinson Islemd Site 
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Fig, 9 Physical horizontal grid point system for 

the free surface model sample problem applied 
to Hutchinson Island site. 
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Dlschargt volume 263,000 G.P.M 
Difcturge Velocity: 0.102 cm/s 
Discharge Temperature: 35.0 C 
Current: 2 cm/ sec 
(Parallel to shoreline) 

Wind: None 

Total Time: 60 min. 



500 m 


Fig. 11 Surface isotherms obtained after 1 hour 
of simulation of the free surface model 
for the sample problem. (Hutchinson Island Site) 
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Air temperature - 20® C 



Figure 12 Wind shear stress reiation. 
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APPENDIX A I 

Wind Stresses I 

I 

The wind shear stresses and are computed by using i 

the Wilson curve as shown in Fig. 12. First, the magnitude 
of the wind velocity, in meters/sec., is used to read off from 
this curve the resultant shear stress^x^^ and x^^ are deter- 
mined by simply resolving x into its respective components. 

As an example, consider a wind of 10 mph from the South 
West direction. Assume that the direction of North is in the 

i 

same direction as the positive x-axis, and East is in the same | 

direction as the positive y-axls as shown below. 



Then t * x cos 45° ■ .4 cos 45° ■ .283 dynes/cm^ 
x^y * X sin 45° ■ .4 sin 45° - .283 dynes/cra"^ 

2 

Where x » .4 dynes/cm for 10 mph » 4.47 meters/sec. 
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APPENDIX B 

HEAT TRANSFER MECHANISMS 

The analysis in this section is taken from Harleman et al. 

( 1973) 

1 . Solar Eradiation (short wave) 

The incident solar radiation impinging on the water surface 
may be expressed as ; 

Where clear sky solar radiation obtained using the 100% 
possible sunshine curve (given in Appendix B) 

C » fraction of sky covered by clouds 
The reflected solar radiation is typically 67o of incident solar 
radiation, hence the net solar radiation absorbed by the water 
surface is: 

•an ■ ”8 ■ ”sr ' 0.94<?g^(l-0.65C^) 

2 . Atmospheric Radiation (long wave) 

The basic equation for the incident atmospheric radiation, 
is given as : 



Where € = average einmitance of the atmosphere 

a = Stefan-Boltzmann constant 
T^* = air temperature (absolute) 
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However, good agreement with experimental data has Indicated 
that Is a function of T. ( ), and specifically, T„*6 

dependence gives best results for atmospheric radiation at low 
temperatures, as well as providing a good fit at high temperatures. 
Clear sky Incident rimospherlc radiation, 9^^ > nay be expressed as: 

and, then Incident atmospheric radiation Including cloudiness 
niay be expressed as : 

H + 0.17C*) 


A figure of 3% Is usually accepted as reflectance of a water 
surface to longwave radiation. Thus the net atmospheric rad- 
iation absorbed by the surface Is : 


P; 


an. 




'ar 


0.97q), 


and, therefore, we have: 


V: 


an 


1.16 X 10“^^CT*)®(1 + 0.17C^) 


3. Longwave Radiation from the Water Surface ^br 

In reference ( ) It Is noted that the emmlsslvlty of a 

water surface Is Independent of temperature and salt or coll- 
oidal concentrations, and gives a value of 0.97. Thus we obtain: 

* 0.97© (T*)"^ 

Where T^ “ water surface temperature. 
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4. Evaporative Heat Flux . ^o 

Evaporation from a water surface occurs as a result of both 
forced (wind driven) convection and free (bouyancy driven) 
convection. The evaporation rom a water surface is usually 
written (mass/area/ time) as: 

E » pr(wj (e. - e.) 

V/here, E « mass flux (mass/area/ time) 

P - density of water 
W » windspeed at height z above surface 
F(Wz)»windspeed function for mass flux including both 
free and forced convection effects (length/ tine/ 
pressure) 

e =» saturated vapor pressure at T 

5 S 

e^“ vapor pressure at height z above surface 
Then writing the above equation in heat units, the evaporative 
heat flux.'^e is given by: 

V/Tiere F(W ) = windspeed function for heat flux (energy/area/ 
tima/pressure) 

Now, dropping the z subscript (and assuming W measured "z" 
above the surface W at the surface) we may expres F(W) for 
a natural water surface and for an artificially heated water 
surface as : 

F(W) = 17W . . . natural water surface 

1/3 

and P (W) =22.4 (Tg-T^) + 14W . . . artificially heated surface 


surface 
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5. Conduction Heat Flux . 

Bowen (see reference ) has suggested that conduction can 
be directly related to evaporative fluxes by assuming that 
eddy dlffuslvltles of heat and mass are Identical. Thus, 


*e " V. 

I 

«h«n Bb ■ 


T -T 

IsJa 


6 -tt 

■ a 


■ Bowen Ratio 


and ■ Bowen constant « 0.255 mm Hg/ F 


and( therefore the conduction heat flux, o , may- be expressed 
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APPENDIX C 

THE EQUILIBRIUM TEMPERATURE AND THE 
SURFACE HEAT TRANSFER COEFFICIENT 
The nec heat transfer ( ) through a water surface is 

composed of radiation penetrating the water surface from above, 
radiation out of the water surface, evaporation, and conduction 
transfer. These are indicated schematically in the following 



Pig. P-1. Heat Transfer Mechanisms at the Water 
Surface . 


The follow^. ig heat balance results, 






'^sr 


^r. " ®e “ 


where 




n 


'sn "^an 

net heat input - cp + ® 

sn 


an 


®br 


- a ...(F-lb) 


N0W| equation (F-1) may be rewritten as, 


<Pn * «Pr 


'Pl 


(P-2) 


Where «p^ * net absorbed radiation * + 9 


an 


and «L “ '^r ®e ®c 


A. 


Equilibrium Temperature Calculation, Te (See Appendix E) 
Under equilibrium conditions equation (F-2) yields. 


*n ' ° ° ’r ■ ®L 
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80 Chat, 

•r “ "t. 


(P-3) 


Then by using the approximate formulae in reference ( ) 

we obtain by setting T. - T_, 

0.94 ’’^(1-0.650*) + X.16 X U+0.17e*) 

■ 0.97,(T,*)® + F(w) C (•,-•,). + Ci,(T,-T^n .. (P-4) 


fiisrc ■ elaar sky solar radiation 

C • cloudiness ratio 
■ air temperature or ®P) 

• equilibrium temperature (®C or ®P) 

^ ■ absolute temperature or ®R) 

F(I7) a> windspeed function (BTU/ftVdayi mm Kg) 
e^ ■* saturated vapor pressure at water surface 
temperature (mm Hg) 

e^^ " saturated vapor pressure at air temperature 
Onm Hg) 

9 m Stefan-Boltzmann constant » 4-1 x .0 BTU/ft , 
day, 

« Bowen constant = 0.255 mm Hg/®P (see Appendix E) 


W • windspeed (mph) 

For a natural water surface, 

F(W) - 17W (F-5a) 

and, for an artificially heated surface, 

F(W) - 22.4 (T -TJ 1/3 + 14W (F-5b) 

6 » 
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Thus, equation (F*4) becomes, 

0.94 9,g(l-0.65C^) + 1.16 X lo"^^ (T^*) ^ (1+0.1 7C^) 
- 0.97<p(T *)^ +• 17W [(•-«.) + 0.255 (T-T)] 

V SO o o 

or , 


0.94 (1-0.650^) + 1.16 x lO*^^ (T^*) ® (1+0.1 7C^) 

■ 0.97(p(T *)^ +[ 22.4 (T-T)^'^^ + 14Vj 

• o o 

t 

•»• 0.255 (T,-T^) ] 

_ - . , cp__i o_, c^, T and VJ •* T can be 

Therefore for known ac a a a. e 

determined by trial and error methods. 

B. Surface Heat Transfer Coefficient (K) 

From reference ( ) the surface heat transfer coefficient 

K, can be determined as follows. 






K • » since « C-(T.) and . . * 0. 

*av av ^ " ®*av 

\dJare *(Tg + T^)/2 


Thua» 


n - 3.88<p(tV + P(W)C(r|^) + C. ] 

av 


+ - ®a> ^ <=b'''s - VI 


av 


^ere 


8F (W) 
3T 

av 



0 

1/3(22.4) 

s a 


for natural water surface 

for artificially heat>^d 
water surface 


1 
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C. Numerical Example 

Consider nacural water surface 

C -0 

T„-25®C 

a 

W >10 mph 

Location - Miami (latitude 26°N) 
Date - December 20 
From reference ( ) • 


•T. - 25®C e. S 0.43 paia 

a a 

•T. ■ 27 C as guess » a 


e ■ 0.51 paia 


a “ — — — . s 

Prom reference ( )# Figure 2,15, pg. 2-61 (soe Figure ?**3) 


i*' a 


^ac 


425 Langleys/day ■ 1560 BTU/ft^/day ••• uaing 

100?teunshine curve at 26^N, Dec. 20 
2 2 

Nbtes 1 L 2 uigley/min. a 220.62 BTU/£t f hr. ■ 1 calarie/cm min. 

« 

Then using equation (F-6) with 0*0, 

0.94 (1560) (1) + 1.16 X 10"^^ (5.37x10^)® (1) » 4250 

4 X lo’® (5.406x10^)^ + 170 C (.255) (2)+(.08) (51.7)] 
4206 close enough! 

.*. «27®C 

(xAxere 1 psia » 51.7 mm Hg) 

Then from equation { ) , 

K -3.88 X 4.1 X 10“® (5.406x10^)^ + 170 (.255+0.0251x51 . 7) 

K 2^290 BTU/ft^, °F, day 


where 


3T 


e -e^ 
^ e a 




= .0251 


T®T 


av 


I 
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D. Ditcuaslon 

The equilibrium surface temperature , T. , for a natural 
water surface, can be greater than the atmospheric temperature. 
T^, whereby T. increases from values below up wO as 
equilibrium is reached. As can be seen in the figure below 
(F*2) T. can be greater or smaller than T depending on the time 
of the day. Simply T.> T. during the hours of sunshine and 

6 A 

> T_ at night when the water surface is cooling. 

A A 



*This plot is taken from ( ) and has no relation to the num 

erical example given in this paper. However . the numerical 
example considered 100% possible hours of sunshine (T > T ) . 

6 A 














